R/km.sim.survtimes.R

In stefkruegel/CPsurv: Nonparametric Change Point Estimation for Survival Data

#' Simulates Survival Times using Kaplan-Meier
#'
#' @inheritParams cpest
#' @param nobs Number of observations.
#' @param weibexp Logical; if \code{TRUE}, survival times above change point have
#'   constant hazard; if \code{FALSE} all survival times are generated by using
#'   the estimated survival curve (relevant for generation of censoring times).
#' @param changeP Change point

km.sim.survtimes <- function(nobs, time, event, weibexp, changeP = NULL){
  survdata <- survival::Surv(time, event, type = "right")

  # Kaplan-Meier estimation
  km <- survival::survfit(survdata ~ 1, type = "kaplan-meier")

  # fitted survival curve
  surv <- summary(km)$surv
  survtimes <- summary(km)$time

  # smoothed survival function
  S <- spline(predict(smooth.spline(round(surv, 4), survtimes),
                      x = seq(1e-04, 1, by = 1e-04)), n=1e4)
  S$x <- round(S$x, 4)
  names(S) <- c("Surv", "survtimes")

  u <- round(runif(nobs, min = 1e-04, max = 0.9999), 4)

  if(weibexp){
    rateE <- sum(event[time > changeP] == 1) / sum(time[time > changeP] - changeP)
    # Survivalfunction at change point
    S.chp <- S$Surv[which.min(abs(S$survtimes - changeP))]

    sim.weib.exp <- function(u, S, S.chp){
      if(u > S.chp) t <- S$survtimes[which(S$Surv == u)]
      else t <- changeP + qexp(p = u/S.chp, rate = rateE, lower.tail = F)
      t
    }

    sapply(u, FUN = sim.weib.exp, S, S.chp)
  } else{
    sapply(u, function(u) S$survtimes[which(S$Surv == u)])
  }
}