I am hoping to simulate survival times from a lung transplant dataset (n > 20,000). In order to do this I need an estimate of the cumulative baseline hazard function.
I have used the
basehaz function in R and then fit a Weibull distribution to it using
nls. However the resulting simulated survival times are significantly lower than the (88% censored) dataset:
The second technique I used is to fit a Weibull distribution (again using
nls) directly to the survival function given by the KM curve. This results in survival curves that match, however I'm not sure of the "correctness" of the survival times since the training dataset is heavily censored and the generated survival times are not censored.
What is the most correct technique to estimate the cumulative baseline hazard function and generate realistic survival times for heavily censored survival data?
Additional details, see response below
I use basehaz with centered = FALSE, would this help address the problem of predictions at the means?
Using survreg I've found that an exponential distribution fits much better than a Weibull distribution, I'm guessing I can't just find which value of lambda matches the survival curve when exp(-lambda * time) is plotted?
I can generate survival times using -(log(U) / (lambda * HR)), where U is a uniform random number in (0, 1] - how would I know if the value of lambda used is appropriate to my survival dataset?
The model has 4 continuous variables, 3 of which are fit with restricted cubic splines and 1 categorical variable with 4 levels.
To simulate the survival times I want to take in a list of patients with the corresponding variables listed above, and then calculate the linear predictor and thus the hazard ratio to generate survival times as in (3) above.
Possible solution, looking for feedback
Thank you very much for the help so far, I feel I have a better idea of what I am doing now. Here is the full process I went through, if it sounds correct please let me know and I can update this post.
I started by using the
cph method from the
rms package and built a Cox model:
rc <- cph(Surv(pdata$waittime, pdata$death_cens) ~ o2_listing + group + rcs(func_status) + rcs(bmi_listing) + rcs(pa_mean), pdata, y=T, x=T)
I then built a dataframe that I planned to use as a reference/baseline:
tmpdf <- data.frame( o2_listing = mean(pdata$o2_listing), group = "A", func_status = median(pdata$func_status), bmi_listing = mean(pdata$bmi_listing), pa_mean = mean(pdata$pa_mean) )
I decided to use the mean value for the continuous variables, the median value for functional status (this is a numeric variable but it takes discrete values) and the reference group "A" for the categorical variable.
I then used
survfit and passed the reference dataframe to it using
newdata and then fit a Weibull distribution to the cumulative hazard function:
sv <- survfit(rc, newdata = tmpdf) DF <- data.frame(time = sv$time, hazard = sv$cumhaz) fit <- nls(DF$hazard ~ l * DF$time ^ v, data=DF, algorithm="port", start=list(l=0.001, v = 1)) l <- coef(fit) v <- coef(fit)
Plotting the survival function for the reference group and overlaying the Weibull curve shows a good fit:
The next step is to calculate the linear predictor (LP) for the reference group:
Survival times are generated using the following function:
(-(log(U))/(l * exp(LP)))^(1/v)
The value of the linear predictor will shift the curve, however since
v were calculated with a reference group with a linear predictor of
-0.4337017 this needs to be subtracted from the LP that is being simulated.
So here is an example, updating the dataframe to specific values:
tmpdf <- data.frame( o2_listing = 6, group = "D", func_status = 50, bmi_listing = 32, pa_mean = 15 )
Again the LP can be calculated:
predict(rc, newdata=tmpdf) which returns
The survival curve can then be plotted using the updated dataframe, then the simulated curve can be overlayed like so:
lines(exp(-l * DF$time^v * exp(0.01615275 - -0.4337017)) ~ DF$time)
This results in the following fit for a new "simulated" individual:
If this appears correct please let me know and I'll update this to hopefully save time for anybody else who needs to simulate survival times.