Following Bender et al. (2003) and
Leemis (1987), simulation of survival
times is possible if there is a function that inverts the cumulative
hazard (\(H^{-1}\)). Random survival
times for a baseline distribution can be generated from a uniform
distribution between 0-1 \(U\) as:
\[ T = H^{-1}(-log(U)) \] For a
survival distribution object, this can be accomplished with the function
rsurv(s_object, n) which will generate n
number of random draws from the distribution s_object. All
objects of the s_distribution family implement a function that inverts
the survival time with the function invCum_Hfx()
The function ggplot_survival_random() helps to graph
Kaplan-Meier graphs and cumulative hazard of simulated times from the
distribution
s_obj <- s_exponential(fail = 0.4, t = 2)
ggplot_survival_random(s_obj, timeto =2, subjects = 1000, nsim= 10, alpha = 0.3)Survival times with hazard proportional to the baseline hazard can be simulated \[ T = H^{-1}\left(\frac{-log(U)}{HR}\right) \] where \(HR\) is a hazard ratio.
The function rsurvhr(s_object, hr) can generate random
numbers with hazards proportional to the baseline hazard. The function
produces as many numbers as the length of the hr vector. for
example:
s_obj <- s_exponential(fail = 0.4, t = 2)
group <- c(rep(0,500), rep(1,500))
hr_vector <- c(rep(1,500),rep(2,500))
times <- rsurvhr(s_obj, hr_vector)
plot(survfit(Surv(times)~group), xlim=c(0,5))
The function
ggplot_survival_hr() can plot simulated data
under proportional hazard assumption.
s_obj <- s_exponential(fail = 0.4, t = 2)
ggplot_survival_hr(s_obj, hr = 2, nsim = 10, subjects = 1000, timeto = 5)Survival times with accelerated failure time to the baseline hazard can be simulated \[ T = \frac{H^{-1}(-log(U))}{AFT}\] where \(AFT\) is an acceleration factor, meaning for example an AFT of 2 has events two times quicker than the baseline
The function rsurvaft(s_object, aft) can generate random
numbers accelerated by an AFT factor. The function produces as many
numbers as the length of the aft vector. for example:
s_obj <- s_lognormal(scale = 2, shape = 0.5)
ggplot_survival_aft(s_obj, aft = 2, nsim = 10, subjects = 1000, timeto = 5)In this example, the scale parameter of the Log-Normal distribution represents the median time, and in this simulation an acceleration factor of 2 moves the median time from 2 to 1
The Proportional Hazards and the Accelerated Failure Time effects can be combined into a single model. Following Chen and Jewell (2001), the Extended Hazards model defines the hazard as \[ h^*(t) = HR \cdot AFT \cdot h_0(AFT \cdot t) \] which gives the cumulative hazard \(H^*(t) = HR \cdot H_0(AFT \cdot t)\) and survival time \[ T = \frac{H_0^{-1}\left(\dfrac{-log(U)}{HR}\right)}{AFT} \] This model nests both models described above as special cases: setting \(AFT = 1\) recovers the Proportional Hazards model, and setting \(HR = 1\) recovers the Accelerated Failure Time model.
Note that this differs from the Accelerated Hazards model of Chen and Wang (2000), \(h^*(t) = h_0(\theta t)\), which rescales the baseline hazard in time without the additional Jacobian factor. That model is obtained from the Extended Hazards model above as the special case \(HR = 1/AFT\) (with \(\theta = AFT\)).
The function rsurveh(s_object, aft, hr) generates random
numbers under the Extended Hazards model. The function produces as many
numbers as the length of the aft/hr vectors,
which must be of the same length. For example:
s_obj <- s_lognormal(scale = 2, shape = 0.5)
ggplot_survival_eh(s_obj, aft = 2, hr = 0.5, nsim = 10, subjects = 1000, timeto = 5)When a subject can have more than one episode over follow-up (e.g. repeated infections or hospitalizations), the time of each subsequent episode can be generated conditional on the time of the previous one, following Leemis (1987). Two assumptions about how risk behaves after an episode give two different processes.
Under a renewal process, risk resets after each episode: the time to
the next episode is a fresh, independent draw from the same (baseline,
proportional hazards, or accelerated failure time) distribution, added
to the time of the previous episode: \[
T_{i+1} = T_i + H^{-1}\left(\frac{-\log(U)}{HR}\right) \] or,
under an accelerated failure time effect, \[
T_{i+1} = T_i + \frac{H^{-1}(-\log(U))}{AFT} \] This is
implemented by renewhr(s_object, hr, prevtime) and
renewaft(s_object, aft, prevtime), which take the time of
the previous episode prevtime and generate the time of the
next one.
Under a non-homogeneous Poisson process, risk does not reset after
each episode: a single cumulative hazard \(H(t)\) accumulates over calendar time from
the start of follow-up, and consecutive episode times satisfy \(H(T_{i+1}) - H(T_i) \sim
\text{Exponential}(1)\), so \[ T_{i+1}
= H^{-1}\left(H(T_i) - \frac{\log(U)}{HR}\right) \] or, under an
accelerated failure time effect, \[ T_{i+1} =
\frac{H^{-1}\left(H(AFT \cdot T_i) - \log(U)\right)}{AFT} \] This
is implemented by nhpphr(s_object, hr, prevtime) and
nhppaft(s_object, aft, prevtime).
The difference between the two processes only shows up when the baseline hazard is not constant. With an increasing baseline hazard, for example, a renewal process keeps generating episodes at the same average pace, since every episode resets the risk back to its starting point. Under a non-homogeneous Poisson process, risk keeps accumulating over calendar time and never resets, so later episodes follow each other increasingly quickly:
# Weibull baseline with an increasing hazard (shape > 1)
s_obj <- s_weibull(scale = 1, shape = 2)
n <- 5000
hr <- rep(1, n)
# First episode, common to both processes
t1 <- rsurvhr(s_obj, hr)
# Renewal process: risk resets at each episode
t2 <- renewhr(s_obj, hr, t1)
t3 <- renewhr(s_obj, hr, t2)
c(gap1 = mean(t1), gap2 = mean(t2 - t1), gap3 = mean(t3 - t2))
#> gap1 gap2 gap3
#> 0.8928635 0.8815392 0.8948339
# Non-homogeneous Poisson process: risk keeps accumulating, never resets
p2 <- nhpphr(s_obj, hr, t1)
p3 <- nhpphr(s_obj, hr, p2)
c(gap1 = mean(t1), gap2 = mean(p2 - t1), gap3 = mean(p3 - p2))
#> gap1 gap2 gap3
#> 0.8928635 0.4388354 0.3335675The average gap between episodes stays roughly constant under the renewal process, while it shrinks with each successive episode under the non-homogeneous Poisson process.