How to simulate survival data with censoring using R

Question

I want to simulate survival data with a sample size of N=100, which follows the Weibull distribution with proportional hazards and constant baseline hazard. Two correlated covariates, which follow different distributions and different censoring rates (e.g. .1, .3, .5) should be integrated. I took the R code from here and varied it a bit, but I only get errors.

R code:

library(survival)
set.seed(123)
simulWeib <- function(N, lambda, rho, beta)
{
correlated covariates X and Y
X <- rexp(N)
  Y <- 2*X + rnorm(N)
Weibull latent event times
v <- runif(n=N)
  Tlat <- (- log(v) / (lambda * exp(X * beta)))^(1 / rho)
censoring times
C <- rnorm(n=20, mean = 0, sd = 1)
follow-up times and event indicators
time <- pmin(Tlat, C)
  status <- as.numeric(Tlat <= C)
data set
data.frame(id=1:N,
             time=time,
             status=status,
             x=X,
             y=Y)
}
fit = coxph(Surv(time, status) ~ x + y, data = dat)

What did I do wrong? A further question is how to divide the censoring into right, left and interval censoring.

This post is related to How to create a toy survival (time to event) data with right censoring.

Answer 1

What did I do wrong?

Your code allows for negative censoring times:

# censoring times
C <- rnorm(n=20, mean = 0, sd = 1)

For a standard Cox model as in your last line of code, negative times aren't allowed. The apparent model for your code sampled from an exponential distribution, with necessarily non-negative values, to get random censoring times. You need to use some distribution similarly limited to non-negative values, for example a log-normal, if you don't want to sample from an exponential distribution..

how to divide the censoring into right, left and interval censoring?

That distinction, in the R survival package, is made in how you set up the Surv() object that represents the survival outcome in the model. See the manual page for Surv(). To avoid confusion when trying to incorporate all types of censoring in a single model, you could use a particular format, Surv(time, time2, event, type="interval2") when submitting your data to your model:

think of each observation as a time interval with (-infinity, t) for left censored, (t, infinity) for right censored, (t,t) for exact and (t1, t2) for an interval. This is the approach used for type = interval2.

If you want to allow for interval censoring you will have to simulate both start and end times for each interval for each individual, making sure not only that both time values are non-negative but also that the end time is greater than the start time for each interval. In practice, you typically have multiple intervals for each individual in that situation, representing for example follow-up times between which an event might have occurred but only detected at time2. Thus for interval censoring you might want to simulate a whole set of follow-up times for each individual and then combine those times in the above format.

You cannot, however, then fit a Cox model to the interval-censored data with the survival package:

Presently, the only methods allowing interval censored data are the parametric models computed by survreg and survival curves computed by survfit; for both of these, the distinction between open and closed intervals is unimportant. The distinction is important for counting process data and the Cox model.

You will have to use another package, like icenReg, for proper Cox modeling of interval-censored data.