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!
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