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

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  • Transformation of the variance-covariance matrix might be possible but it seems likely to create many opportunities for error. The matrix produced originally is (asymptotically) multivariate normal; after transformation it won't be. In general, it's best and simplest to do all calculations in the original scale and only at the end to transform to the desired scale for display. As an example, in a Cox model you work with the coefficients and their covariances in the original log-hazard scale and only at the end do the exponentiation to the hazard-ratio scale. Why complicate matters? – EdM May 4 at 13:16 • So for the time being, do you think it's reasonable to take the default log-linear scale parameter from survfit() via coef(wFit), transform it to the original scale of Weibull via scale <- exp(coef(wFit)), and then take the variance-covariance matrix from survfit() via vcov(wFit) as-is without transformation? I need the var-cov matrix to run MASS::mvrnorm(). May 4 at 15:05 • No. What I have in mind is, until the very very end, to work completely with the variance-covariance matrix from the standard survreg() or corresponding flexsurvreg() fit and the corresponding coefficient estimates (via model$icoef for survreg() or coef(model) for flexsurvreg()). That's the scale in which the coefficient estimates are asymptotically multivariate normal so that mvrnorm() can be used reliably. The survfit() function doesn't fit models (except Kaplan-Meier curves); it generates survival curves from a previously fitted Cox or accelerated-failure-time model.
– EdM
May 4 at 15:19
• Sorry for the sloppy question, everywhere I used survfit() I meant to use survreg(). I'll try your recommendation. May 4 at 15:33

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))^exp(-survregCoefs))
}

# 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)