1

I am using R's flexsurvreg function (in the flexsurv package) to fit a AFT model to my data.

This is the line of code that fits the model to the data:

TestModel <- flexsurvreg(Surv(time,death) ~ param1 + param2 + param3 + param4 + param5 + param6 + param7 + param8 + param9 + param10 + param11 + param12 + param13, data = DataTest, dist = "weibull")

Once the model fits, this is a summary of the results:

Estimates: 
        data mean     est        L95%       U95%       se         exp(est)   L95%       U95%     
shape      NA         9.99e-01         NA         NA         NA         NA         NA         NA
scale      NA         2.20e+02         NA         NA         NA         NA         NA         NA
param1     1.32e-01   2.51e-01         NA         NA         NA   1.29e+00         NA         NA
param2     1.61e-01  -1.54e-02         NA         NA         NA   9.85e-01         NA         NA
param3     1.89e-01  -4.68e-02         NA         NA         NA   9.54e-01         NA         NA
param4     1.76e-01  -2.25e-02         NA         NA         NA   9.78e-01         NA         NA
param5     1.87e-01  -5.35e-02         NA         NA         NA   9.48e-01         NA         NA
param6     7.56e-01  -2.74e-01         NA         NA         NA   7.60e-01         NA         NA
param7     2.28e-01   3.23e-02         NA         NA         NA   1.03e+00         NA         NA
param8     1.58e-01  -1.69e-02         NA         NA         NA   9.83e-01         NA         NA
param9     4.32e-01  -1.89e-02         NA         NA         NA   9.81e-01         NA         NA
param10    1.30e+02  -1.01e-03         NA         NA         NA   9.99e-01         NA         NA
param11    2.26e+01  -4.08e-03         NA         NA         NA   9.96e-01         NA         NA
param12    5.54e+02  -2.84e-04         NA         NA         NA   1.00e+00         NA         NA
param13    9.57e+01  -4.69e-03         NA         NA         NA   9.95e-01         NA         NA

N = 40320,  Events: 32154,  Censored: 8166
Total time at risk: 2584693
Log-likelihood = -171611.5, df = 15
AIC = 343253.1

I want to measure how the covariates affect the survival time. The estimates provide an understanding of this. Also, as I read here, $exp(est)$ provides an estimate of how the hazard changes with change in 1 unit of a covariate by keeping the other covariates fixed. Is there a way I can calculate p-values for these covariates?

I have fitted a Weibull distribution to my dataset.

Community
  • 1
statBeginner
  • 1,251
  • 2
  • 17
  • 22
  • 1
    Can you add the line code you're using with the `flexsurvreg()` function? At the moment the output (from `summary` presumably on the fitted object) doesn't have any standard errors printed, which seems odd, and these would need to be calculated before any statistical inference could be done (confidence intervals, hypothesis tests) – James Stanley May 14 '15 at 22:22
  • I have modified the post. Also, the distribution is Weibull, not exponential as in the original post. – statBeginner May 14 '15 at 22:53
  • 2
    The issue with not having standard errors and confidence intervals in the summary output still stands; without them, you can't do the statistical inference you want. Try first doing a standard coxph rather than the flexsurvreg; also note that the Weibull shape parameter estimate is essentially 1, so you may just have an exponential distribution. I suspect there's a problem with presenting your data correctly to flexsurvreg, although I don't have experience with that function myself. – EdM May 15 '15 at 14:16

1 Answers1

3

There are no intervals and standard errors printed, and I assume that is because the model fit did not converge, therefore the reported estimates are not the real maximum likelihood estimates. There should have been a message from flexsurvreg saying this. Non-convergence can happen when you are trying to fit a model with too many parameters in it to a dataset that is too small, so a solution is to use fewer covariates.

If the same model worked in "survreg" from the "survival" package, then tweaking the optimisation options used by flexsurvreg may help. Adding

flexsurvreg(..., control=list(ndeps=rep(1e-06, 15)))

(replacing 15 if necessary by the number of parameters in your model) worked for a user with a similar problem. flexsurvreg appeared to be having difficulty getting the Hessian, and consequently the SEs and CIs.

Also a point of misunderstanding, it's not a Cox model, it's a Weibull parametric model, and it's an accelerated failure time model, not proportional hazards. See the vignette in the package for further explanation.

Chris Jackson
  • 176
  • 3
  • I am sorry. Yes, it is an accelerated failure model. My bad. Here is the problem: Since flexsurvreg internally calls survreg for fitting the parametric weibull, the odd thing I notice is survreg does provide me with the standard errors for the same data and same covariates. – statBeginner May 19 '15 at 21:28
  • 2
    This shouldn't happen. As the author of flexsurvreg, could you send me the data (anonymised if necessary) and code to reproduce this, and I'll investigate. – Chris Jackson May 19 '15 at 22:03
  • That is a problem Chris. The data is confidential. I will get in touch with the authorities and see how it can be managed and get back to you soon. – statBeginner May 20 '15 at 19:09
  • Could you try adding an argument control = list(ndeps=1e-06, 15) to flexsurvreg? This tightens the step size used to approximate the Hessian. This worked for a user with a similar problem, where survreg worked but flexsurvreg didn't. flexsurvreg appeared to be having difficulty getting the Hessian, and consequently the SEs and CIs. – Chris Jackson Jun 02 '15 at 09:31
  • When I add the argument, i get an error that states: Error in optim(control = list(ndeps = 1e-06, 15), method = "BFGS", par = c(-0.000836815874312455, : 'ndeps' is of the wrong length In addition: Warning message: In optim(control = list(ndeps = 1e-06, 15), method = "BFGS", par = c(-0.000836815874312455, : unknown names in control: – statBeginner Jun 10 '15 at 04:36
  • 1
    Sorry, I meant list(ndeps = rep(1e-06, 15)) – Chris Jackson Jun 12 '15 at 15:11
  • Right. That solves the problem. Can you include that as part of the answer please? – statBeginner Jun 12 '15 at 16:35