Generating Random Survival Times From Any Hazard Function

Written by Peter Rosenmai on 14 Apr 2017.

Let's get 1,000 random survival times (for use, perhaps, in a simulation) from a constant hazard function (hazard = 0.001):

hazard_fn = function(t) rep(0.001, length(t))
survival_times = random_survival_times(hazard_fn, 1000)

And let's check that the Kaplan-Meier curve for these survival times appproximates, as expected, the curve P(t) = exp(-0.001t):

require("survival")
plot(survfit(Surv(survival_times, rep(1, length(survival_times)))~1), xlab="Time", ylab="Survival Probability",
     main="Sampled and Expected Survival Curves for h(t) = 0.001")
neg_exp_fn = function(x){exp(-0.001 * x)}
curve(expr=neg_exp_fn, from=0, to=max(survival_times), add=TRUE, col="red", lwd=2)

Kaplan-Meier curve for 1,000 survival times sampled from h(t) = 0.001

Good. That looks much as we'd expect.

This works with any hazard function. So here's another example, a hazard function that has hazard rise linearly with age:

hazard_fn = function(t) 0.001 + t * 0.0001
survival_times = random_survival_times(hazard_fn, 1000)
require("survival")
plot(survfit(Surv(survival_times, rep(1, length(survival_times)))~1), xlab="Time", ylab="Survival Probability",
     main="Sampled Survival Curve for h(t) = 0.001 + t * 0.0001")

Kaplan-Meier curve for 1,000 survival times sampled from h(t) = 0.001 + t * 0.0001

Finally, here's the code that does the work:

inverse = function(fn, min_x, max_x){
   # Returns the inverse of a function for a given range.
   # E.g. inverse(sin, 0, pi/2)(sin(pi/4)) equals pi/4 because 0 <= pi/4 <= pi/2
   fn_inv = function(y){
      uniroot((function(x){fn(x) - y}), lower=min_x, upper=max_x)[1]$root
   }
   return(Vectorize(fn_inv))
}

integrate_from_0 = function(fn, t){
   int_fn = function(t) integrate(fn, 0, t)
   result = sapply(t, int_fn)
   value  = unlist(result["value",])
   msg    = unlist(result["message",])
   value[which(msg != "OK")] = NA
   return(value)
}

random_survival_times = function(hazard_fn, n, max_time=10000){
   # Given a hazard function, returns n random time-to-event observations.
   cumulative_density_fn = function(t) 1 - exp(-integrate_from_0(hazard_fn, t))
   inverse_cumulative_density_fn = inverse(cumulative_density_fn, 0, max_time)
   return(inverse_cumulative_density_fn(runif(n)))
}