3

I am trying to estimate the mean time to failure for a Weibull distribution fitted to some survival data with flexsurvreg from the flexsurv package. I need to be able to estimate the standard error for use in a simulation model.

Using flexsurvreg with the lung data as an example;

require(flexsurv)
lungS <- Surv(lung$time,lung$status)
lungfit <- flexsurvreg(lungS~1,dist="weibull")
lungfit

Call:
flexsurvreg(formula = lungS ~ 1, dist = "weibull")

Maximum likelihood estimates: 
est   L95%   U95%
shape   1.32   1.14   1.52
scale 418.00 372.00 469.00

N = 228,  Events: 165,  Censored: 63
Total time at risk: 69593
Log-likelihood = -1153.851, df = 2
AIC = 2311.702

Now, calculating the mean is just a case of plugging in the estimated parameter values into the standard formula, but is there an easy way of getting out the standard error of this estimate? Can survreg do this?

B--rian
  • 5,578
  • 10
  • 38
  • 89
peedeerich
  • 253
  • 3
  • 10

1 Answers1

4

In flexsurv version 0.2, if x is the fitted model object, then x$cov is the covariance matrix of the parameter estimates, with positive parameters on the log scale. You could then use the asymptotic normal property of maximum likelihood estimators. Simulate a large number of multivariate normal vectors, with the estimates as means, and this covariance matrix (using e.g. rmvnorm from the mvtnorm package). This gives you replicates of the parameter estimates under sampling uncertainty. Calculate the corresponding mean survival for each replicate, then take the SD or quantiles of the resulting sample to get the standard error or a confidence interval.

Chris Jackson
  • 871
  • 5
  • 7
  • Thanks Chris, I eventually figured out that this approach would work well. Remember if you are doing this to use the logged parameter mean vector, i.e. `x$res.t` not `x$res` (not that I would make such an elementary mistake - ha!). I guess I was being lazy in hoping for a canned solution, but the simulation took much less time than asking the question! – peedeerich Aug 15 '13 at 20:07
  • Just a note that version 0.5 of the package has some new utilities for doing this kind of thing. You can produce summaries for user-defined functions of the parameters by passing the "fn" argument to "summary.flexsurvreg", which uses the (visible) function "normboot.flexsurvreg" to generate a sample from the asymptotic normal distribution of the estimates. – Chris Jackson Oct 24 '14 at 09:16