In the code shown at the bottom of this post, I plot survival curves for the lung
dataset from the survival
package using a fitted exponential model, using the K-M nonparametric model, and run/show simulations using the exponential model.
I use bootstrapping, resampling from the original data with replacement to create multiple bootstrap samples using sample()
. For each bootstrap sample, the code fits the exponential distribution using the survreg()
function. This process is repeated, generating a distribution of estimates, representing the variability and uncertainty of the exponential statistical model.
My objective with this ultimately is given a partial survival curve (say 500 periods of the lung dataset), generating conservative simulations for periods 501-1000. I don't show that in this code example. When drafting similar code for the Weibull distribution, I use both bootstrapping (with sample()
function) and additionally simulated uncertainty of the Weibull parameters using MASS:mvrnorm()
, to derive a nicely dispersed range of simulation outcomes.
However, in this exponential model example, the exponential distribution has only one parameter, the rate (λ) parameter; so MASS:mvrnorm()
makes no sense in this case. To introduce more dispersion in outcomes in the below code I use rnorm(1, mean = 0, sd = 0.05)
in the sim_params
section (all commented out in the code and in the below illustration to not introduce this additional uncertainty factor), which as the code is currently drafted is subjective (by manually inputting the SDEV value) and not grounded in the actual data unlike my use of MASS:mvrnorm()
for the Weibull distribution.
So my questions are (1) is there a way to ground this parameter uncertainty factor (sim_params...
) in the actual lung
data? and (2) is this method of modeling uncertainty both using bootstrapping with sample()
and modeling uncertainty in the distribution parameters themselves (in the sim_params
section) theoretically valid?
The image below only shows the results of running the code with only bootstrap resampling functioning, and showing a run of 2000 simulations:
Code:
library(survival)
num_simulations <- 2000
# Fit the exponential model to the dataset
fit <- survreg(Surv(time, status) ~ 1, data = lung, dist = "exponential")
time <- seq(0, 1000, by = 1)
# Compute the exponential survival function using fitted model
survival <- 1 - pexp(time, rate = 1 / exp(fit$coef))
# Generate bootstrap samples and fit exponential models to each sample
bootstrap_fits <- lapply(1:num_simulations, function(i) {
sample_data <- lung[sample(nrow(lung), replace = TRUE), ]
fit <- survreg(Surv(time, status) ~ 1, data = sample_data, dist = "exponential")
return(fit)
})
# Generate random distribution parameter estimates for simulations
sim_params <- sapply(bootstrap_fits, function(fit) {
rate <- fit$coef
params <- rate # this is a bypass of "perturbation" below
# perturbation <- rnorm(1, mean = 0, sd = 0.05) # Adjust sd for simulation dispersion
# perturbed_rate <- rate + perturbation
# params <- perturbed_rate
return(params)
})
# Compute the survival curves for each simulation using the sampled parameters
sim_curves <- sapply(
1:num_simulations,
function(i) 1 - pexp(time, rate = 1 / exp(sim_params[i]))
)
plot(time, survival, type = "n", xlab = "Time", ylab = "Survival Probability",
main = "Survival Plot of Lung Dataset")
sim_lines <- data.frame(
time = time,
do.call(cbind, lapply(1:num_simulations, function(i) {
curve <- sim_curves[, i]
lines(time, curve, col = "lightblue", lty = "solid", lwd = 0.25)
return(curve)
})))
colnames(sim_lines)[-1] <- paste0("surv", 1:num_simulations)
# Compute and add to the plot the Kaplan-Meier survival curve for the dataset
lines(survfit(Surv(time, status) ~ 1, data = lung), col = "blue", lwd = 1)
# Plot the exponential survival curve
lines(time, survival, type = "l", xlab = "Time", ylab = "Survival Probability", col = "red", lwd = 3)
legend("topright",
legend = c("Fitted exponential model",
"Kaplan-Meier & confidence intervals",
"Simulations"),
col = c("red", "blue", "lightblue"),
lwd = c(3, 1, 0.25),
lty = c(1, 1, 1), # 1 = solid, 2 = dashed
bty = "n")