# How to simulate survival data with censoring using R

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.

• Where does your Y go? Commented May 5, 2021 at 14:29
• @Jarle Tufto Y is defined as Y <- 2*X + rnorm(N) or what do you mean in detail?
– Tino
Commented May 5, 2021 at 15:19
• What do you want to accomplish when you “divide the censoring into right, left and interval censoring“? What type of situation are you trying to simulate that way?
– EdM
Commented May 6, 2021 at 3:06
• @EdM I want to simulate all three situations of censoring (right, left and interval) in one model.
– Tino
Commented May 6, 2021 at 7:34

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.

• Thank you for your good answer! I have changed the censoring time back to an exponential distribution. But where can I change the censoring rate (relative share of censored individuals)? What means rexp(n=20) for the censoring time in detail? Is the Y missing in the definition of Tlat <- (- log(v) / (lambda * exp(X * beta)))^(1 / rho)? If yes, where should I add it?
– Tino
Commented May 6, 2021 at 13:04
• @Tino the formula for censoring time in the model code is C <- rexp(n=N, rate=rateC), where N is the total number of cases and rateC determines how quickly censoring drops off with time. So your formula's use of rexp(n=20) is also incorrect here. Changing rateC effectively changes the fraction of cases censored, as a case is censored if C<T.
– EdM
Commented May 6, 2021 at 13:58
• @Tino The model you used for simulation assumes that only X is associated with outcome so that Y is a correlated predictor that has no independent association with outcome. If you want Y also to be associated with outcome on its own, you need to specify separate $\beta$ values for each of X and Y in your simulation function and write $\exp (\beta_x X +\beta_Y Y)$ instead of $\exp(X * beta)$ in the formula for Tlat.
– EdM
Commented May 6, 2021 at 14:06