I am interested in learning about how to simulate survival times from a Survival Model (e.g. Cox-PH).

For example, suppose we fit a Cox-PH Model on some data in R:

# Fit the Cox model
fit <- coxph(Surv(time, status) ~ age + sex + ph.ecog, data = lung)

The results of the model look something like this:

coxph(formula = Surv(time, status) ~ age + sex + ph.ecog, data = lung)

             coef exp(coef)  se(coef)      z        p
age      0.011067  1.011128  0.009267  1.194 0.232416
sex     -0.552612  0.575445  0.167739 -3.294 0.000986
ph.ecog  0.463728  1.589991  0.113577  4.083 4.45e-05

Likelihood ratio test=30.5  on 3 df, p=1.083e-06
n= 227, number of events= 164 
   (1 observation deleted due to missingness)

I have seen R packages that are able to perform similar tasks (e.g. https://www.jstatsoft.org/article/view/v097i03) - but is it possible to write a step-by-step process in R that can simulate from such a model?

I would be interested in learning how to do this.



2 Answers 2


Understanding the Cox PH Model

The survival function under the proportional hazards assumption is

$$ S(t \mid X) = S_0(t)^{\exp(X\beta)} $$

So there are two parts you need to specify:

  • The log hazard ratios $\beta$, and
  • The baseline or reference survival function $S_0(t)$.

Once you have these, you can leverage the fact that $S(t\mid X) = 1-F(t\mid X)$ to draw survival times. Here, $F$ is the CDF of failure times.

Let's walk through this step by step.

Specifying the log hazard ratios

This is a fairly simple step. I'm going to simulate 2 continuous variables and 1 binary indicator. The continuous variables will have mean 0 and standard deviation 1. The binary indicator will have a prevalence of 0.5. I'll specify the effects of the continuous variables to be 1 and 0.325 on the log hazard ratio scale, and the effect of the binary variable to be 0.15 on the log hazard ratio scale.


# No intercept in the cox model
N <- 25000
age <- rnorm(N)
sex <- rbinom(N, 1, 0.5)
weight <- rnorm(N)
X <- model.matrix(~age + sex + weight -1)
beta <- c(1, 0.15, .325)
xb <- X %*% beta

Specifying $S_0(t)$

The cox model is a semi-paremetric model, so it doesn't really matter what we choose here. I'm going to use a weibull distribution so that

$$S_0(t) = \begin{cases}e^{-(x / \lambda)^k}, & x \geq 0 \\ 0, & x<0\end{cases}$$

for appropriate choice of scale and shape. This will allow us to compare predicted survival curves against the true survival curves.

Drawing Survival Times

This is the tricky bit. $S_0(t)$ is the survival, which means $F(t) =1-S_0(t)$ is the CDF. To sample the survival times, we can

  • Draw a uniform random variable, $u$. If we wanted to draw from the PDF of survival times, we would solve the equation $F(t) =u$ for t. This can be done by inverting the CDF or by root finding..
  • Well, if $S_0(t)^{\exp(X\beta)}$ is our survival function, then $F(t \mid X) = 1-S_0(t)^{\exp(X\beta)} $ is our CDF, so we can solve $F(t \mid X) = u$ to get a draw from the survival time distribution.

We can also add censoring if we like. An event is considered censored if the failure time is larger than the censoring time. For simplicity, let's draw censoring times from the same distribution as the survival times.

finv = \(x, u, xb) (1-pweibull(x, 4, 5))^(exp(xb)) - u

failure_times <- rep(0, N)
censor_times <- rep(0, N)

for(i in 1:N){
  failure_times[i] <- uniroot(finv, interval=c(0, 20), u=uu, xb=xb[i])$root
  censor_times[i] <- uniroot(finv, interval=c(0, 20), u=cc, xb=xb[i])$root


time <- pmin(failure_times, censor_times)
event <- as.integer(failure_times<censor_times)

Fit the Model

Now, we just do the model fitting and compare against the parameters we specified

fit <- coxph(Surv(time, event) ~ age + sex + weight)

coxph(formula = Surv(time, event) ~ age + sex + weight)

  n= 25000, number of events= 12464 

           coef exp(coef) se(coef)      z Pr(>|z|)    
age    0.986221  2.681085 0.011754 83.909   <2e-16 ***
sex    0.151580  1.163671 0.017949  8.445   <2e-16 ***
weight 0.319033  1.375797 0.009291 34.337   <2e-16 ***
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

       exp(coef) exp(-coef) lower .95 upper .95
age        2.681     0.3730     2.620     2.744
sex        1.164     0.8593     1.123     1.205
weight     1.376     0.7269     1.351     1.401

Concordance= 0.73  (se = 0.002 )
Likelihood ratio test= 7969  on 3 df,   p=<2e-16
Wald test            = 7474  on 3 df,   p=<2e-16
Score (logrank) test = 7533  on 3 df,   p=<2e-16

The immense sample size if the simulation ensures our estimates are very close to the truth, which they are.

Plotting $S_0(t)$ and The Known Survival Function

Like I said, $S_0$ corresponds to subjects who have 0 for all their covariates. We've cooked up this example so that $1-S_0(t)$ is the weibull CDF. We can plot the estimated survival curve against the weibull survival curve to prove this to ourselves.

fit <- coxph(Surv(time, event) ~ age + sex + weight)

sfit = survfit(fit, newdata = list(age=0, sex=0, weight=0))

  estimated = sfit$surv,
  truth = 1 - pweibull(sfit$time, 4, 5)
) %>% 
  gather(which, surv, -time) %>% 
  ggplot(aes(time, surv, color = which)) + 
  geom_step(size=2, alpha = 0.5)

The two lines are nearly ontop of one another enter image description here

Edit: How Good Was Our Estimate

We can bascially treat our estimated survival curve as a prediction and compute MSE and RSME

  estimated = sfit$surv,
  truth = 1 - pweibull(sfit$time, 4, 5)) %>% 
    mse = mean((estimated - truth)^2),
    rmse = sqrt(mean((estimated - truth)^2))

# A tibble: 1 × 2
         mse    rmse
       <dbl>   <dbl>
1 0.00000900 0.00300

Or compute the pointwise relative error

  estimated = sfit$surv,
  truth = 1 - pweibull(sfit$time, 4, 5)) %>% 
    rel_err = abs(truth-estimated)/truth
  ) %>% 
  ggplot(aes(time, rel_err)) + 

Bigger means worse approximation. I'm not surprised the tail is so bad, the probabilities are very small so any small change results in a large relative

enter image description here

  • $\begingroup$ @ Demetri Pananos: Thank you so much for your answer! Just a question - you assumed that the Event Times have a Weibull Distribution. You then simulated new times from a Weibull Distribution and fit a new model to these simulated times. You then compared the results of the original model vs the new model and concluded that both models were similar. Is there a way to quantify this comparison? Thank you so much! $\endgroup$
    – stats_noob
    Feb 7 at 4:11
  • $\begingroup$ You can look at the pointwise relative error between the estimated survival function and the true survival function. Alternatively, use any of the accuracy metrics (e.g. RMSE) from machine learning $\endgroup$ Feb 7 at 4:21
  • $\begingroup$ @ Demetri Pananos : Thank you for your reply! An idea I had: I was looking at the last picture you created which shows the original survival curve and the new survival curve literally on top of each other. Perhaps there might be some statistical test (e.g. log rank) that can compare how similar the original survival curve and the new survival curve are? $\endgroup$
    – stats_noob
    Feb 7 at 4:23
  • $\begingroup$ @stats_noob I don't see a benefit of a statistical test in this particular case. What is the null in this case: that the two are the same? We know the pink curve estimates the blue curve, there is no to test that. $\endgroup$ Feb 7 at 4:28
  • 1
    $\begingroup$ @stats_noob Done. Just treat the estimate like a prediction from an ML algorithm. $\endgroup$ Feb 7 at 4:41

This answer shows the general principle, and illustrates with a Weibull baseline hazard that follows proportional hazards. You sample randomly from a uniform distribution on $(0,1)$, and find the time that corresponds to that quantile of the survival curve, given the covariate values. That works in general.

Under proportional hazards, the survival over time as a function of time-constant covariate values in a vector $X$ and corresponding regression coefficients $\beta$ is:

$$S(t|X)= S_0(t)^{\exp(\beta' X)}, $$

where $S_0(t)$ is the baseline survival function. So, given your randomly sampled survival fraction and your covariates $X$, you invert the survival function to get the corresponding survival time.

That poses problems with an empirical Cox model. The baseline survival function is a step function, constant between event times and taking sharp steps down at event times. So it's a good idea at least to smooth the empirical baseline survival function. Furthermore, survival data often have a right-censored value for the last observation time, so in a semi-parametric Cox model there is no information on survival beyond that last observation time. That's why you will more typically see simulations based on parametric survival functions.

  • $\begingroup$ @ EdM: thank you so much for your answer! $\endgroup$
    – stats_noob
    Feb 7 at 4:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.