3
$\begingroup$

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.

enter image description here

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)
  }
)
$\endgroup$
1
  • 1
    $\begingroup$ This appears to be a very roundabout way to go about it. And you are assuming that the bootstrap distribution mimics the sampling distribution which is often not the case. $\endgroup$ Jun 2, 2023 at 13:38

1 Answer 1

2
$\begingroup$

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.

$\endgroup$
4
  • $\begingroup$ By paucity I mean working with data when part way through a study. In the lung example the death rate over the complete study time window (max time 1022) is 72.4% calculated as the number of status = 2 out of 228 patients. If we were sitting in the 500th month of the study, 60.5% of the patients would have died. Based on the data through month 500, how can we start molding our expectation of future death rate probabilistically? Going back to earlier queries like stats.stackexchange.com/questions/614198/… $\endgroup$ Jun 4, 2023 at 5:53
  • $\begingroup$ I note your answer in stats.stackexchange.com/questions/608652/… where I could use binomial for live/die, and then estimate the time curves for deaths through survival or other. However, I prefer survival because in addition I'll need to introduce multi-state variables such as stage 1, 2, 3, 4, as patients get sicker which I understand the survival package accommodates. $\endgroup$ Jun 4, 2023 at 5:58
  • 1
    $\begingroup$ @Village.Idyot unless you have complete follow up on all patients, those calculations of "death rates" are erroneous. The denominator continues to include people lost to follow up for reasons other than death. If you do have complete follow up that might be OK, but that's not typical of survival studies and certainly is not the case in the lung data set. Survival analysis takes that problem into account, provided that censoring isn't informative about death status. If your real data have complete follow up, why not ask a question directly about analyzing your real data? $\endgroup$
    – EdM
    Jun 4, 2023 at 6:16
  • 1
    $\begingroup$ @Village.Idyot for example, if you believe the lognormal model, the expected fraction of individuals who will have died by time=500 is 69.2%, and by time=1022 it's 87.6%. Those are both substantially higher than your estimates for those times based on observed deaths and the original numbers of patients. You are correct that more events leads to better precision in fitting a model, but the real problem in predicting beyond the last observation time is whether you have fit the proper parametric form. Perfectly identifying coefficients of the wrong model accomplishes nothing for prediction. $\endgroup$
    – EdM
    Jun 4, 2023 at 13:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.