Motivation

In many scenarios, we may be interested in simulating survival times from a population based on the corresponding life tables. For instance, in relative survival analysis (Rubio et al. 2019) (Eletti et al. 2022) (Basak et al. 2026) simulation studies require the simulation of times to event based on life tables. These life tables are stratified on certain characteristics such as age, year, and deprivation level. The life table defines a piecewise-constant rate hazard function, for each combination of characteristics, which can be used to simulate the times to event.

The SimLT R package allows for simulating times to event, based on the information in a life table.

The following R code presents a step-by-step illustrative example om the simulation of times to event based on life tables from England for the years 2010-2015. The SimLT R package requires functions from the R packages msm and lubridate, which I decided not to include as dependencies. Thus, the user needs to install them and load them manually.

See also: GHSurv, LBANS HazReg

Illustrative example: R code

In this example we will simulate \(n=1000\) times to event using a life table from England based on year, sex, deprivation level, age, English region (gor). The life table can be downloaded from SimLT.

Data preparation

rm(list=ls())

# Required packages
#pak::pak("FJRubio67/SimLT")
library(SimLT)
library(msm)
library(lubridate)
## 
## Attaching package: 'lubridate'
## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union
#------------------------------------------------------------------------
# Life tables
#------------------------------------------------------------------------
LT <- as.data.frame(read.table("England_LT_2010_2015_dep_gor.txt",header = T))
head(LT)
##   X_year sex dep age      rate gor
## 1   2010   1   9   0 0.0047365  10
## 2   2011   1   9   0 0.0046188  10
## 3   2012   1   9   0 0.0046188  10
## 4   2013   1   9   0 0.0046188  10
## 5   2014   1   9   0 0.0046188  10
## 6   2015   1   9   0 0.0046188  10

Simulated design matrix

We now simulate the individuals’ characteristics using the built-in function simDesMatrix. These characteristics will be used to simulate the times to event. This function is based on real data on colon and lung (male and female) patients in England. You can use your own design matrix instead.

#------------------------------------------------------------------------
# Design matrix simulation 
#------------------------------------------------------------------------
# Sample size
n <- 1000
# Design matrix
Xcf <- simDesMatrix(seed = 123, n = n, admin.cens = "31-12-2015", scale.age = FALSE,
                    site = "colon", sex = "female")
head(Xcf)
##   ID  date.diag dep gor      age fup.time
## 1  1 2010-05-11   2   2 80.09770 5.639973
## 2  2 2010-12-21   2   3 90.00683 5.026694
## 3  3 2010-03-31   1   2 76.93449 5.752225
## 4  4 2010-03-11   3   6 92.35372 5.806982
## 5  5 2010-12-17   3   6 67.86141 5.037645
## 6  6 2010-05-13   4   7 28.95104 5.634497

Follow-up period

Define the follow-up period. This requires commands from the R package lubridate.

#------------------------------------------------------------------------
# Follow-up period
#------------------------------------------------------------------------
# fixed points (new year and last day of follow-up)
year.init <- 2010
last.day <- as.Date("2015-12-31")
fixed.times <- c(seq(as.Date("2010-01-01"), as.Date("2015-01-01"), by = "year"),last.day)

Time points and hazard rates

Now, we need to define the time points that define the step function associated to the piecewise-constant rate hazard function.

#------------------------------------------------------------------------
# Time points and hazard rates
#------------------------------------------------------------------------

# Date of diagnosis and date of birth (29 February is moved to 1 March)
dod <- as.Date(Xcf$date.diag)
dob <- as.POSIXlt(dod - Xcf$age*365.25)
leap <- dob$mday == 29 & dob$mon == 1
dob$mday[leap] <- 1
dob$mon[leap] <- 2
dob <- as.Date(format(dob, "%Y-%m-%d"))

# sex
sex0 <- 2

# Candidate time points for all patients at once:
# date of diagnosis, new years, birthdays, and last day of follow-up
yrs <- year.init:year(last.day)
cand <- data.frame(
  id   = c(1:n, rep(1:n, each = length(yrs)), rep(1:n, each = length(yrs)), 1:n),
  date = c(dod,
           rep(as.Date(paste0(yrs, "-01-01")), n),
           as.Date(paste0(rep(yrs, n), rep(format(dob, "-%m-%d"), each = length(yrs)))),
           rep(last.day, n)))
# Keep the points between diagnosis and the last day of follow-up, sorted by patient and date
cand <- cand[cand$date >= dod[cand$id] & cand$date <= last.day, ]
cand <- unique(cand)
cand <- cand[order(cand$id, cand$date), ]

# Age (completed years) at each time point
cand$age <- year(cand$date) - year(dob[cand$id]) -
  (format(cand$date, "%m%d") < format(dob[cand$id], "%m%d"))

# Extracting the hazard rates from the life table for each patient at "age + t":
# a single hashed lookup on (year, sex, dep, age, gor) for all patients
LT.key <- paste(LT$X_year, LT$sex, LT$dep, LT$age, LT$gor)
key <- paste(year(cand$date), sex0, Xcf$dep[cand$id], cand$age, Xcf$gor[cand$id])
cand$rate <- LT$rate[match(key, LT.key)]

# Dates, lengths of times in days between time points, and hazard rates for each patient
dates  <- unname(split(cand$date, cand$id))
times  <- lapply(dates, function(d) c(0, as.numeric(diff(d))))
hrates <- unname(split(cand$rate, cand$id))

# List containing the time points, hazard rates, and dates
Xcf_rates <- list(times = times, hrates = hrates, dates = dates)

Simulating the times to event

Now, we simulate the times to event using the command sim_pophaz. This requires a seed for the simulation, and a list with time points at which the piecewise hazard jumps (times) and the values of the piecewise hazard (hrates). This list was constructed in the previous subsection.

##################################################################################################
# Population survival times
##################################################################################################
simC <- sim_pophaz(seed = 123,lst = Xcf_rates, method = "direct")
sim.pop <- simC$sim.pop
status.pop <- simC$status

# Censoring
table(status.pop)
## status.pop
##   0   1 
## 794 206
library(survival)

# Kaplan-Meier estimator for the simulated times
kmsim <- survfit(Surv(sim.pop, status.pop) ~ 1)

# Comparison
plot(kmsim$time, kmsim$surv, type = "l", col = "black", lwd = 2, lty = 1, ylim = c(0,1),
     xlab = "Time", ylab = "Survival", cex.axis = 1.5, cex.lab = 1.5, main = "")

Basak, P., C. Maringe, F. J. Rubio, and A. R. Linero. 2026. “Understanding Inequalities in Cancer Survival Using Bayesian Machine Learning.” Journal of the American Statistical Association: Applications and Case Studies 121 (553): 72–84.
Eletti, A., G. Marra, M. Quaresma, R. Radice, and F. J. Rubio. 2022. “A Unifying Framework for Flexible Excess Hazard Modelling with Applications in Cancer Epidemiology.” Journal of the Royal Statistical Society: Series C (Applied Statistics) 71: 1044–62.
Rubio, F. J., L. Remontet, N. P. Jewell, and A. Belot. 2019. “On a General Structure for Hazard-Based Regression Models: An Application to Population-Based Cancer Research.” Statistical Methods in Medical Research 28: 2404–17.