2
$\begingroup$

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.

Improve this question
$\endgroup$
3
  • 1
    $\begingroup$ 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) $\endgroup$
    – James Stanley
    May 14, 2015 at 22:22
  • $\begingroup$ I have modified the post. Also, the distribution is Weibull, not exponential as in the original post. $\endgroup$
    – statBeginner
    May 14, 2015 at 22:53
  • 2
    $\begingroup$ 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. $\endgroup$
    – EdM
    May 15, 2015 at 14:16

1 Answer 1

Reset to default
3
$\begingroup$

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.

Improve this answer
$\endgroup$
7
  • $\begingroup$ 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. $\endgroup$
    – statBeginner
    May 19, 2015 at 21:28
  • 2
    $\begingroup$ 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. $\endgroup$
    – Chris Jackson
    May 19, 2015 at 22:03
  • $\begingroup$ 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. $\endgroup$
    – statBeginner
    May 20, 2015 at 19:09
  • $\begingroup$ 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. $\endgroup$
    – Chris Jackson
    Jun 2, 2015 at 9:31
  • 1
    $\begingroup$ Sorry, I meant list(ndeps = rep(1e-06, 15)) $\endgroup$
    – Chris Jackson
    Jun 12, 2015 at 15:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.