# Simulate time to event times based on an existing subset of data

I have to task to predict data volume to be processed for a study that currently has 40 patients, but will have 100 patients eventually. Since most of the data is generated during the treatment phase, i would like to focus on this phase.

For these 40 patients, i calculated the time to end of treatment (censoring patients still on treatment). I thought that using this survival curve, i would be able to 'draw' another 60 patients from the same curve. If i were to have the survival times (e.g. time to end of treatment) for all 100 patients I would be able to add up all time to end of treatment times, to get an idea of the total number of months of treatment in our sample. To do this, i would assume for example 5 patients to be accrued per month. Then add the simulated survival time for each patient, take a close out date and sum 'treatment time' for all patients.

I have googled everywhere and found all kinds of packages that can simulate the most complex survival data, however, i fail to find/understand how to simulate 'simple' survival data.

Does anyone know how i can simulate X survival times using existing data? Surely i can do it roughly by hand as there are many knowns (e.g. estimated survival time for various proportions of the data) but surely there would be an easier method?

The lung dataset might provide a good example...

 sfit <-survfit(Surv(time,status)~1,data=lung)
plot(sfit)


This dataset contains 228 patients. How would you simulate data for 1000 patients with similar characterisitics of the first 228 patients (e.g. median survival time of 310 days, etc)?

EDIT: i guess i just need to be able to sample from the Weibull (or exponential) distributions. Question is, how to estimate the Weibull parameters using the data that i already have, in a way that it incorporates cencored data....

• There's no one right way to simulate anything. You really haven't asked a specific programming question here (just because you want to do the simulation in R, doesn't automatically make that a programming question). If you are unfamiliar with methods that will simulate data with the properties you desire, you should probably post to Cross Validated for statistical advice instead. – MrFlick Sep 23 '14 at 4:13
• Should have been posted on CrossValidated. The request to use "existing data" is unclear unless it is a request to deploy bootstrap methods. The rms/Hmisc combo package has power analysis strategies and worked examples of simulation on the cph and cpower help pages. – DWin Sep 23 '14 at 6:18
• In the Cran Task View on Survival Analysis, there is an entire section on Simulation... Perhaps you find something there. – shadow Sep 23 '14 at 7:26
• Anyone??? I guess I just need to be able to sample from the Weibull (or exponential) distributions. Question is, how to estimate the Weibull parameters using the data that i already have (or from the Kaplan-Meier curve), in a way that it incorporates cencored data.... – Luc Sep 25 '14 at 1:57

## 1 Answer

You should first decide on a survival time distribution that best fits the data you have already, this can be done by fitting several parametric distributions (e.g exponential, Weibull, log-normal e.t.c the R package 'flexsurv' will be useful as it provides the AIC as part of its return values) to the data from the 40 patients and comparing the AIC of the models.

Below is a mock of how to choose between say a log-normal and an exponential model for the lung cancer data.

library(flexsurv)
surv<-with(lung, Surv(time, status))  # create the survival object
model<-flexsurvreg(surv ~ 1, dist="lnorm")  # fit the log-normal model
model
Call:
flexsurvreg(formula = surv ~ 1, dist = "lnorm")
Estimates:
est    L95%   U95%
meanlog  5.660  5.510  5.820
sdlog    1.100  0.983  1.230

N = 228,  Events: 165,  Censored: 63
Total time at risk: 69593
Log-likelihood = -1169.269, df = 2
AIC = 2342.538


Lets have a visual of the fit

plot(model, ylab="Survival probability", xlab="Time")
legend("topright",legend=c("KM Plot","Fitted"), lty=c(1,1),col=c("black","red"), cex=0.5)


Now lets try the Weibull model

model2<-flexsurvreg(surv ~ 1, dist="weibull") # fit weibull model
model2
Call:
flexsurvreg(formula = surv ~ 1, dist = "weibull")
Estimates:
est     L95%    U95%
shape    1.32    1.17    1.49
scale  418.00  372.00  469.00

N = 228,  Events: 165,  Censored: 63
Total time at risk: 69593
Log-likelihood = -1153.851, df = 2
AIC = 2311.702


See how it fits below

plot(model2, ylab="Survival probability", xlab="Time")
legend("topright",legend=c("KM Plot","Fitted"), lty=c(1,1),col=c("black","red"), cex=0.5)


Model 2 (the Weibull distribution) fits better, it has a lower AIC. This is also supported by the graph which shows smaller deviations between the fitted and observed failure probabilities over time for the Weibull distribution. You can then generate your random samples from your choice of the best fitting distribution. I have used Weibull but it might not be the case with your data.

n<-100 # number of samples to draw
T<-rweibull(n, shape=1/exp(model2$coefficients['shape']), scale = exp(model2$coefficients['scale'])) #  time to event
C<-rweibull(n, shape=1/exp(model2\$coefficients['shape']), scale = 450)  # random time to decide censoring, scale parameter above that estimated from the sample data
(sum(T>C)/n)*100 # censoring rate (%), adjust the scale parameter for 'C' to make this as close as possible to the observed censoring rate in the sample data

T_new<-NULL
for (i in 1:n) T_new[i] <- min(T[i],C[i]) # the final time to use as your sampled time
d <- (T<C)*1  # censoring indicator: event occurs if d==1


Now you have the time and the censoring indicator, you actually don't need the later from your description; your interest is in the total of T_new. You will actually have to repeat picking samples of 100 several times, say 1000, then pick the average of the total of T_new across the replications.

I hope this was helpful.