Does this approach to simulation for survival analysis, of breaking the analysis into deaths versus survivors, appear reasonable?

Question

I've spent last several weeks learning about survival analysis, see one of the last posts at How to simulate variability (errors) in fitting a gamma model to survival data by using a generalized minimum extreme value distribution in R?

Now I am primarily concerned with simulating death rates and secondarily deriving survival curves for the deceased. Ultimately, this is leading towards simulating deaths/survival using extreme value distributions with heavy right tails (even if not best-fitting) for simulating conservative, very-bad-case scenarios especially when dealing with a paucity of data. The code below is a first step in that direction.

Does the approach I describe and per the code below appear reasonable? If so, are there easier or better approaches?

1. I use the lung dataset from the survival package as my example.
2. I use bootstrap sampling (bootSample() in code below) to derive death rates (deathRate <- ...) and to extract the lung data for only the deaths from the same bootstrapped samples where "status" == 2 (bootDeaths[[i]] <<- ...).
3. Using AIC, lognormal provided the best fit for bootstrap sampled death rates. Code not shown for this goodness-of-fit testing.
4. I draw lognormal random samples for each of the bootstrap samples and derive a histogram of death rates per the image below on the left.
5. I then take the deaths from the same bootstrapped samples and fitting the deaths with the survreg() function and applying the lognormal distribution, plot their survival curves (plot_survival_curves(...)) as shown in the image below on the right.

Code:

library(MASS)
library(survival)

nbr <- 100
timeLine <- seq(0, max(lung$time))
bootDeaths <- list()

# Use bootstrapping for both average death rates and for plotting survival curves for deaths
bootSample <- sapply(
  1:100,
  function(i) {
    sampleData <- lung[sample(nrow(lung), replace = TRUE), ]
    bootDeaths[[i]] <<- sampleData[sampleData$status == 2, ] # used in plotting death survival curves later
    deathRate <- with(sampleData, mean(status == 2))
    return(deathRate)
  }
)

### Generate random samples for the lognormal distribution, calculate and plot death rates ###
fit <- MASS::fitdistr(bootSample,"lognormal")
params <- fit$estimate
sampLognorm <- rlnorm(1000, params[1], params[2])
hist(sampLognorm, breaks = "FD", col = "steelblue",
     xlab = "Death rate", ylab = "Frequency", main = "Histogram of Lognormal Samples")
sampDeathRate <- mean(bootSample)
abline(v = sampDeathRate, col = "black", lty = 1, lwd = 3)
popDeathRate <- with(lung, mean(status == 2))
abline(v = popDeathRate, col = "red", lty = 1, lwd = 3)
legend("topright", legend = c(paste("Sample Average:", round(sampDeathRate, 4)),
                              paste("Population Average:", round(popDeathRate, 4))),
       lty = c(1,1), lwd = c(3,3), col = c("black", "red"), bty = "n")

### Lognormal survival curves for patients who die ###
plot(timeLine, type = "n", xlab = "Time", ylab = "Survival Probability", main = "Lung Data Survival Plot", ylim = c(0, 1), xlim = c(0,max(lung$time)))

# Fit lognormal distribution and plot survival curves for each deceased sample
plot_survival_curves <-  sapply(
  1:nbr,
  function(i){
    sampleDat <- data.frame(bootDeaths[[i]])
    fit <- survreg(Surv(time, status == 2) ~ 1, data = sampleDat, dist = "lognormal")
    meanlog <- fit$coef  
    sdlog <- fit$scale  
    surv_prob <- 1 - plnorm(timeLine, meanlog = meanlog, sdlog = sdlog)
    lines(seq(0,length(surv_prob)-1), surv_prob, col = "lightblue", lty = "solid", lwd = 0.25)
  }
)

Answer 1

One big problem comes to mind: even without censored observations, there is no single "death rate" (e.g., deaths per person at risk per year) unless there's an exponential survival curve. The hazard function is the continuous death rate over time, and it's far from constant over time for a lognormal distribution; see the NIST page for example plots of lognormal hazard functions.

Without an exponential survival curve, any "death rate" based on dividing a number of deaths by a number at risk will depend on the particular time window involved. What exactly do you mean by a "death rate" in this context? How would you apply it in practice?

In your situation this problem is exacerbated by omitting right-censored observations from the calculations. That necessarily introduces bias into the estimate, in a way that depends heavily on the censoring pattern in the underlying data. It's hard to think of a scenario in which omitting censored event times leads to anything other than trouble.

Survival analysis done properly allows for other than exponential survival curves while handling censored event times. It allows for predictions of survival at specific times of interest, so that if you are particularly concerned about, say, early events you can focus attention on them.

In terms of the sampling variability you want to capture (variability in survival-curve estimates arising from random sampling of the population) the covariance matrix of the coefficient estimates should contain the information you need.

If the data set is large enough and the form of the model is adequate, then the asymptotic multivariate normality of the estimates can be used directly without resampling. I showed that at the end of my answer to the question you linked, where resampling from the survival distribution gave essentially the same coefficient covariance as the original model's coefficient covariance matrix.

If the form of the model isn't adequate, then modeling on resampled cases might give an estimate of a type of variability: the variability in coefficient estimates in a model that doesn't properly fit the data. What's the point? You are better served by developing an adequate model.

If the data set isn't large enough, you have to consider whether to trust the model at all. Your statement about "paucity of data" is thus troubling. It might be more useful to present the actual scenario that you are interested in, as I fear that the "paucity of data" will limit what can be accomplished.