extract_raw_sim.R

#' Functions to extract raw simulated samples
#'
#' @name extractrawsim
NULL

#' @rdname extractrawsim
#' @export
#' @param sim A `survparamsim` class object generated by [surv_param_sim()] function.
#' @details
#' [extract_sim()] extracts raw survival time & event status for all simulated subjects.
extract_sim <- function(sim) {

  time.var   <- as.character(attributes(stats::formula(sim$survreg))$variables[[2]][[2]])
  status.var <- as.character(attributes(stats::formula(sim$survreg))$variables[[2]][[3]])
  if(status.var[[1]] == "!") status.var <- status.var[[2]]


  if(methods::is(sim, "survparamsim_resample")){
    sim.merged.with.cov <-
      sim$newdata.nona.sim %>%
      dplyr::select(-dplyr::all_of(time.var), -dplyr::all_of(status.var), -n.resample) %>%
      dplyr::left_join(sim$sim, ., by = c("rep", "subj.sim")) %>%
      dplyr::select(rep, subj.sim, time, event, dplyr::everything())

  } else if(methods::is(sim, "survparamsim_pre_resampled")){
    sim.merged.with.cov <-
      sim$newdata.nona.sim %>%
      dplyr::select(-dplyr::all_of(time.var), -dplyr::all_of(status.var)) %>%
      dplyr::left_join(sim$sim, ., by = c("rep", "subj.sim")) %>%
      dplyr::select(rep, subj.sim, time, event, dplyr::everything())

  } else {
    sim.merged.with.cov <-
      sim$newdata.nona.sim %>%
      dplyr::select(-dplyr::all_of(time.var), -dplyr::all_of(status.var)) %>%
      dplyr::left_join(sim$sim, ., by = c("subj.sim")) %>%
      dplyr::select(rep, subj.sim, time, event, dplyr::everything())
  }

  return(sim.merged.with.cov)
}


#' @rdname extractrawsim
#' @export
#' @param hr.pi a return object from [calc_hr_pi()] function.
#' @details
#' [extract_hr()] extracts simulated HR for all repeated simulations.
#' If HR was calculated based on simulated survival times with
#' [calc_hr_pi()] function, it also returns p values for Cox regression
#' fits, one for each group
#' based on Wald test and another for the overall significance of the
#' coefficient based on logrank test. The latter has the same values
#' across treatment groups when >2 levels in treatment
extract_hr <- function(hr.pi) {

  return(dplyr::select(hr.pi$sim.hr, -description))
}


#' @rdname extractrawsim
#' @export
#' @param km.pi A return object from [calc_km_pi()] function.
#' @details
#' [extract_km_obs()] extracts observed Kaplan-Meier curves.
extract_km_obs <- function(km.pi) {
  return(km.pi$obs.km)
}

#' @rdname extractrawsim
#' @export
#' @param km.pi A return object from [calc_km_pi()] function.
#' @details
#' [extract_medsurv()] extracts simulated median survival times for all repeated simulations
extract_medsurv <- function(km.pi) {
  if(methods::is(km.pi, "survparamsim.kmpi.aveHR")){
    stop("Median survival calculation not implemented in `calc_ave_km_pi()`` yet")
  }
  return(km.pi$sim.median.time)
}

#' @rdname extractrawsim
#' @export
#' @param km.pi A return object from [calc_km_pi()] function.
#' @details
#' [extract_medsurv_delta()] extracts delta of median survival times between treatment groups
extract_medsurv_delta <- function(km.pi) {

  pi.range   <- km.pi$pi.range
  trt.assign <- km.pi$trt.assign
  sim.median.time <- km.pi$sim.median.time

  if(length(km.pi$trt.syms) == 0) stop("`trt` needs to be specified in `calc_km_pi()`")
  if(length(km.pi$trt.syms) > 1) stop("`trt` can only take one string in `calc_km_pi()")

  group.syms <- km.pi$group.syms
  trt.sym    <- km.pi$trt.syms[[1]]

  # Check trt values
  check_trt(km.pi$median.pi, trt.sym)

  # Convert trt to factor
  sim.median.time <-
    sim.median.time %>%
    dplyr::mutate(!!trt.sym := factor(!!trt.sym))

  # Reverse control vs trt
  if(trt.assign == "reverse"){
    sim.median.time <-
      sim.median.time %>%
      dplyr::mutate(!!trt.sym := forcats::fct_rev(!!trt.sym))
  }

  sim.median.time <-
    sim.median.time %>%
    dplyr::mutate(.trt.group.index = as.integer(!!trt.sym),
                  .trt.group.index = paste0(".trt.group.", .trt.group.index))

  trt.group.index.map <-
    sim.median.time %>%
    dplyr::select(!!trt.sym, .trt.group.index) %>%
    dplyr::distinct()

  # Calculate delta
  sim.median.time.delta <-
    sim.median.time %>%
    # Change to wide data
    dplyr::select(!(!!trt.sym)) %>%
    dplyr::select(!n) %>%
    tidyr::pivot_wider(names_from = .trt.group.index,
                       values_from = median) %>%
    # Calc delta for all the non-control groups
    dplyr::rename(.trt.control.group = .trt.group.1) %>%
    dplyr::mutate(dplyr::across(dplyr::starts_with(".trt.group."), ~.x-.trt.control.group)) %>%
    dplyr::select(!.trt.control.group) %>%
    # Convert back to long data and recover the original trt variables
    tidyr::pivot_longer(dplyr::starts_with(".trt.group."),
                        names_to = ".trt.group.index",
                        values_to = "median_delta") %>%
    dplyr::left_join(trt.group.index.map, by = ".trt.group.index") %>%
    dplyr::select(!.trt.group.index)

  # Reverse back control vs trt
  if(trt.assign == "reverse"){
    sim.median.time.delta <-
      sim.median.time.delta %>%
      dplyr::mutate(!!trt.sym := forcats::fct_rev(!!trt.sym))
  }

  sim.median.time.delta %>%
    dplyr::select(rep, !!!group.syms, !!trt.sym, dplyr::everything()) %>%
    dplyr::arrange(rep, !!!group.syms, !!trt.sym) %>%
    return()
  }