Options for transforming the variance-covariance matrix generated by the survreg() function to the original scale of the Weibull distribution?

Question

I'm working with the survreg() function of the R survival package, and I understand that the default scale parameter for the Weibull distribution generated by this function is on the log-linear scale. In the code snippet at the bottom of this post, you can see I transform the scale parameter to the original scale of the Weibull distribution with exp(coef(wFit)). I also need to extract the variance-covariance matrix for this Weibull fit, and my understanding is the variance-covariance matrix would also need be transformed to the original scale of the Weibull distribution. From my research I have found different methods for transforming the variance-covariance matrix and I am confused as to whether these methods are correct, and which I should apply for my needs.

Basically, I am trying to forecast survival probabilities into future periods. For testing purposes, I truncate a survival curve (as you can see in post How to generate multiple forecast simulation paths for survival analysis? where the lung1 object is a hypothetical truncation to 500 periods of the lung dataset) and then plot a curve (or curves) for future periods via the fit of the partial curve. Conservative estimates are better; better to understate survival probabilities than it is to overstate.

Please, any guidance on the below attempts, listed as Method A, B, and C? Are any of them correct or advisable, given my objectives described above? Using the lung dataset from survival as the basis for experimentation. Pointing me in the direction of any digestible reference materials will also help!

Another method that I have seen is simulation, involving simulating many datasets based on the estimated coefficients and scale parameters, transforming them to the original scale of the Weibull distribution, and then calculating the variance-covariance matrix of the transformed parameters. I'm not ready to try simulation until I better understand the simpler methods!

Code:

library(survival)

wFit <- survreg(Surv(time,status)~1, dist="w", data=lung)
scale <- exp(coef(wFit))
shape <- 1/wFit$scale

### Method A ###
transform_mat1 <- vcov(wFit) * (scale^2)
transform_mat1

### Method B ###
# Outline the transformation matrix
transform_mat2 <- matrix(c(1, 0, 0, -1/scale^2), nrow = 2)

# Apply transformation to the new transform_mat matrix
cov_est_transformed <- transform_mat2 %*% vcov(wFit) %*% t(transform_mat2)
cov_est_transformed

### Method C ###
# Extract the log-linear scale parameter estimates from the model
log_estimates <- coef(wFit)
log_scale <- log_estimates  # the second parameter is log(scale)

# Compute the gradient of the log-scale parameter estimate with respect to the original-scale parameters
d_log_scale_d_scale <- 1 / exp(log_scale)  # partial derivative of log(scale) w.r.t scale

Answer 1

Below is how I take EdM's advice per his comments above, where I work with the variance-covariance matrix using the standard survreg() format and transform to the original scale using the weibCurve() function per the below code. See for complete explanation: How to generate multiple forecast simulation paths for survival analysis?

Code:

library(survival)

simNbr <- 10

# converts from survreg() log-linear scale to standard parameterization used by dweibull()
weibCurve <- function(time, survregCoefs) {
  exp(-(time / exp(survregCoefs[1]))^exp(-survregCoefs[2]))
}

# Fit the Weibull model to the dataset
fit <- survreg(Surv(time, status) ~ 1, data = lung, dist = "weibull")

# Generate bootstrap samples and fit Weibull models to each sample
bootstrap_fits <- lapply(
  1:simNbr, 
  function(i) {
    sample_data <- lung[sample(nrow(lung), replace = TRUE), ]
    fit <- survreg(Surv(time, status) ~ 1, data = sample_data, dist = "weibull")
    return(fit)
  }
)

# Generate random Weibull parameter estimates for simulations
simParams <- sapply(bootstrap_fits, function(fit) {
  newCoef <- fit$icoef
  params <- MASS::mvrnorm(1, mu = newCoef, Sigma = vcov(fit))
  return(params)
})

# Simulate uncertainty in the parameters of the Weibull distribution and plot results
simPaths <- sapply(
  1:simNbr, 
  function(i) {
    params <- simParams[, i]
    return(weibCurve(time, params))
  }
)

# Create the plotSims dataframe and label column headers
plotSims <- data.frame(time = time, simPaths)
colnames(plotSims)[-1] <- paste0("surv", 1:simNbr)