Probabilistic estimates of Kaplan-Meier curves using S(t) and std.err for each time point

Question

I am working on a survival analysis where I am fitting different parametric models to survival data. Varying the models that I have fit to the data is straightforward, via Cholesky decompositions, but I'd like to generate probabilistic Kaplan-Meier curves using the "survival" and "std.err" parameters shown in the data below (I am working in R).

Call: survfit(formula = Surv(futime, fustat) ~ 1, data = PF)

time n.risk n.event survival std.err lower 95% CI upper 95% CI
    0    221       0    1.000  0.0000       1.0000        1.000
    1    149      52    0.764  0.0286       0.7104        0.823
    2     50      24    0.587  0.0396       0.5147        0.670
    3     33       3    0.541  0.0445       0.4608        0.636
    4     26       1    0.522  0.0469       0.4377        0.623
    5     16       4    0.428  0.0577       0.3288        0.558
    6      9       3    0.336  0.0653       0.2300        0.492
    7      8       1    0.299  0.0679       0.1917        0.467
    8      6       2    0.224  0.0685       0.1233        0.408
    9      2       2    0.140  0.0642       0.0571        0.344
   10      2       0    0.140  0.0642       0.0571        0.344
   11      2       0    0.140  0.0642       0.0571        0.344

One idea that I've had is to use the estimated survival probability and variance for each time point to sample from a normal distribution for each time point. There is an obvious problem with this – I could end up sampling a higher survival probability for time 3 than for the time before it, which is illogical.

If I enforce crude "rules" (i.e. force the next sampled $S(t)$ to be equal to or less than the previous $S(t)$ estimate) I will end up systematically underestimating $S(t)$.

I feel as though I am barking up the wrong tree and I'm pretty stumped. Has anyone come across a solution for this before? I'd be grateful for any insights.

Answer 1

One trick would be to sample from the distribution of the hazard function $h(t)$ itself rather than from the integrated hazard function that underlies the survival function $S(t)$. For a Kaplan-Meier analysis this would only be done at discrete observed event times. At event time $i$ the estimated hazard $\hat h_i$ is $d_i/n_i$, where $d_i$ is the number of events at that time and $n_i$ is the number still at risk. The variance of $\hat h_i$ at time $i$ is then $h_i (1-h_i)/n_i$. (See the derivation of Statistical considerations on the Wikipedia Kaplan-Meier page.) You would start for example with the 221 cases initially at risk and adjust both the numbers of events and the numbers at risk at each event time accordingly (although I haven't thought enough about this to say how to deal properly with censoring in this situation).

If you are doing parametric modeling it would make more sense to sample from the underlying event probability distribution. For example, a Weibull model has a corresponding extreme value distribution of log-survival times:

If the Weibull has the shape parameter γ and the characteristic life α, then the extreme value distribution (after taking natural logarithms) has μ = ln(α), β=1/γ.

So sampling from the corresponding extreme value function would be a reasonable approach, less tied to the vagaries of your initial data set. Other parametric survival models have corresponding probability distributions for survival time, as explained for example in these course notes.