Correct technique to estimate cumulative baseline hazard function in R?

Question

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

1. I use basehaz with centered = FALSE, would this help address the problem of predictions at the means?

2. 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?

3. 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)[1]
v <- coef(fit)[2]

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:

predict(rc, newdata=tmpdf)

-0.4337017

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 l and 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 0.01615275.

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.

Answer 1

The problem is just what the "baseline" for the "baseline hazard function" is. If you define that carefully, you will be OK. The way that you describe what you've done suggests that you fell into a trap. From the help page for basehaz():

This function is simply an alias for survfit, which does the actual work and has a richer set of options. The alias exists only because some users look for predicted survival estimates under this name...If there are factor variables in the model, then the default predictions at the "mean" are meaningless since they do not correspond to any possible subject; correct results require use of the newdata argument of survfit.

Thus I suspect that what basehaz() returned didn't represent any realistic combination of covariate values. As the manual page suggests, it's safest to specify your own set of covariate values via the newdata argument to survfit() so that you know exactly what your baseline represents.

Furthermore, you can't safely compare the raw Kaplan-Meier curve to any specific baseline hazard, regardless of which covariate values you specify for the baseline. The mean of survivals (one crude way to think about the Kaplan-Meier curve) is not the same as survival at mean covariate values, given the inherent nonlinearity of survival models. That also means there is a danger in using a Weibull baseline estimate based on the raw Kaplan-Meier curve.

If you want to work with a Weibull baseline hazard, why not just fit a full Weibull survival model directly? The R survreg()function does that quite simply. It will return the location and scale parameters for a defined "baseline" model and the coefficients that represent the change in location as covariate values change.

In response to edited question

If you specify centered = FALSE to basehaz() then you get the baseline hazard with all covariate values set to 0. That can make it easier to set up calculations for any particular set of covariate values than you have with the default of centered = TRUE. The behavior of basehaz() can depend on the specific software used; a quick test suggests that cph and coxph objects return different results with the default. It might be that cph objects internally specify centered = FALSE, but you should check the documentation.

For parametric modeling, it's not clear why you "found that an exponential distribution fits much better than a Weibull distribution." The exponential is a specific case of the Weibull, with the scale parameter (in the location-scale format) forced to take a value of 1. I wonder if there's some other problem here. An exponential or Weibull model provides an (Intercept) that corresponds to baseline values of the covariates, which you then can translate to a rate constant. With the R survreg() function and an exponential model, the exponential rate constant $\lambda = \text{exp}(-(\text{Intercept})) $. Be careful, particularly with Weibull, as parameterizations vary.

If you have the correct baseline value of $\lambda$ at reference covariate values, then your formula for simulation seems OK. You can save yourself a lot of trouble by using the R simsurv package, which is designed for this purpose. For censored data you also have to simulate censoring times for all cases, then choose the lower of the event or censoring time as the time outcome, labeled accordingly.