Periodic jumps in piece-wise constant hazard simulation

[jgboh]

Hello,

when I use a piece-wise constant hazard function to simulate data I always get periodic jumps after the first cutpoint, resulting from accumulated eventtimes. In the first image below this is most visible in the 4th intervall. The survial curve should look like resulting from an exponential distribution since the hazard is constant.

I used my own code as well as the example in section 4.4. provided in the package authors publication When I reduce the number of cut-points that separate the hazard function to 4 the simulated data shows these jumps.

library(simsurv)
library(survival)

##### create a bathtub-shaped hazard function
set.seed(1729)
# ncuts <- 19
ncuts <- 4 
cuts <- sort(rexp(ncuts, rate = 0.1)) 

pw_times <- c(0, cuts)

N <- length(pw_times)
pw_haz <- sort(abs(rnorm(N)))

pw_haz <- abs(pw_haz - median(pw_haz))

##### user defined hazard function 
haz <- function(t, x, betas, lb_interval, haz_interval, ...) {
  haz_interval[findInterval(t, lb_interval)]
}

##### Simulate data
cov <- data.frame(id = 1:5000)
dat <- simsurv(x = cov, hazard = haz, lb_interval = pw_times,
               haz_interval = pw_haz)
summary(dat$eventtime)

##### Plot the piecewise constant hazard function
plot(x = pw_times, y = pw_haz, type = "s",main = "Piecewise constant hazard function", xlab = "Time", ylab = "Hazard rate")

##### Plot KM curve of simulated data
fit_km <- survfit(Surv(eventtime, status) ~ 1 ,data = dat)
plot(fit_km, conf.int = FALSE, main = "Kaplan-Meier curve of simsurv simulated data", xlab = "Time", ylab = "S(t)")
abline(v = cuts, col = "grey", lty = 2)

I run into the same issue withmy own piecewiese constant hazard function. here there are is only one cutpoint at t=5. In the second intervall a treatment effect is added that seperates both arms.

### piecewise constant hazard function 
pc_hazard <- function(t, x, betas) {
  ### calculate hazard at time t
  ifelse(test = (t <= 5), 
         yes = betas[["beta1"]], 
         no = betas[["beta2"]] * x[["treatment"]])
}
### Set simulation parameters
set.seed(1234)
n <- 1000
t_max = 10
# Covariate data
covs <- data.frame(id = 1:n,
                   treatment = rep(c(1,2), each = (n/2)))

# Population (fixed effect) parameters
betas <- c(0.02, 0.1)
names(betas) <- c("beta1", "beta2")

#---------------------------------------------------------------------------------
### Simulation
sim_times <- simsurv(hazard = pc_hazard, 
                     x = covs, 
                     betas = betas, 
                     maxt = t_max)

sim_data <- merge(sim_times, covs)

#---------------------------------------------------------------------------------
### Visualisation of the simulated data
fit_km <- survfit(Surv(eventtime, status) ~ treatment, data = sim_data)
plot(fit_km,     
     conf.int = FALSE, 
     col = c("red", "blue"), 
     main = "Simulated delayed effect", xlab = "t", ylab = "S(t)")