sim_pexp_cr: Simulate competing risks time-to-event data via piece-wise...
In adibender/pem.xgb: GBT (PEM)

Piece-wise exponential implementation of simulation algorithm described in Beyersmann et al. 2009 (<doi: 10.1002/sim.3516>).

1	sim_pexp_cr(formula, data, cut)

`formula`	An extended formula that specifies the linear predictor. If you want to include a smooth baseline or time-varying effects, use `t` within your formula as if it was a covariate in the data, although it is not and should not be included in the `data` provided to `sim_pexp`. See examples below.
`data`	A data set with variables specified in `formula`.
`cut`	A sequence of time-points starting with 0.

library(pammtools)
library(dplyr)
library(survival)
# Example from Competing Risks book for comparison
simul.dat.cp <- function(n, h01, h02, cens.param) {
  times <- rexp(n, h01 + h02)
  ev <- rbinom(n, size = 1, prob = h01 / (h01 + h02))
  ev <- ifelse(ev == 0, 2, 1)
  cens.time <- runif(n, cens.param[1], cens.param[2])
  obs.time <- pmin(times, cens.time)
  obs.cause <- as.numeric(times == obs.time) * ev
  data.frame(obs.time, obs.cause)
}
set.seed(923)
n <- 1000
dat.exo7 <- simul.dat.cp(n, 0.5, 0.9, c(0.5, 3))
# We compute the Nelson-Aalen estimates using survfit()
# That’s just to get the number of events and the risk set
temp01 <- survfit(Surv(obs.time, obs.cause == 1) ~ 1, dat.exo7)
temp02 <- survfit(Surv(obs.time, obs.cause == 2) ~ 1, dat.exo7)
na01 <- cumsum(temp01$n.event / temp01$n.risk)
na02 <- cumsum(temp02$n.event / temp02$n.risk)

plot(temp01$time, na02, type="s", ylab="Cumulative transition hazard", xlab="time")
lines(temp01$time, na01, type="s", col=2)

# create data set with variables which will affect the hazard rate
# (not used here, but usually more complex examples of interest)
df <- cbind.data.frame(x1 = runif (n, -3, 3), x2 = runif (n, 0, 6)) %>%
 as_tibble()
set.seed(24032018)
df <- cbind.data.frame(
  x1 = runif (n, -3, 3),
  x2 = runif (n, 0, 6))
# two component formula specifying cause specific hazards
form <- ~ log(0.5)| log(0.9)
sim_df <- sim_pexp_cr(form, df, seq(0, 3, by =.25)) %>%
 mutate(
  cens_time = runif(n(), 0.5, 3),
  status = if_else(cens_time < time, 0, 1),
  time = pmin(time, cens_time),
  type = status * type)
temp01_2 <- survfit(Surv(time,type == 1) ~ 1, sim_df)
temp02_2 <- survfit(Surv(time, type == 2) ~ 1, sim_df)
na01_2 <- cumsum(temp01_2$n.event / temp01_2$n.risk)
na02_2 <- cumsum(temp02_2$n.event / temp02_2$n.risk)

lines(temp01_2$time, na02_2, type="s", col = 3)
lines(temp01_2$time, na01_2, type="s", col = 4)