Closed Geoff-Holmes closed 6 years ago
Hi,
This looks weird and I've not experienced problems like this with the Weibull distribution. (In a previous version of survHE
there was an error in the code for the Gompertz distribution, which would make the calculation of the survival curves all wrong --- but this has been long fixed). In all my tests with both parameterisations of the Weibull, things have worked out just fine...
Can you share the data so I can reproduce the error and debug?
Unfortunately, I can't share the data that I used for this. However, as soon as I can I will see if I can adapt it or reproduce the problem in some way that I can share.
Thanks for the email correspondence on this Gianluca.
So, if I understand correctly, this is down to the way survHE parses the formula passed to it.
I noticed previously that whilst in flexsurv, I can do:
survObj<-with(dat, Surv(Survival, event))
formula<-survObj~group
m<-flexsurvreg(formula, data=dat, dist="weibull")
In survHE, this doesn't work (Error in attributes(terms(formula))$variables[[2]][1] : object of type 'symbol' is not subsettable
) and the call to Surv has to be included explicitly in the passed formula.
This is not a problem, but I think coming from familiarity with flexsurv, it tripped me up somewhat.
In the current issue, the formula I specified was:
formula<-with(dat, Surv(Survival, Status==1)~Group)
This formula works okay in flexsurv, but for survHE, the event needs to first defined as a binary:
dat$Event<-dat$Status==1
And then,
formula<-with(dat, Surv(Survival, Event)~Group)
m2<-fit.models(formula, data=regdata, distr="weibull", method="hmc")
plot(m2)
all works fine and agrees with m1
in the original post.
Yes. In fact, you can check this by comparing the two version of the model (MLE or HMC) when using a binary indicator Event
(taking value 1 if fully observed and 0 if right censored --- irrespective of how the right censoring comes about) with the version of the MLE (flexsurv-based) model when you simply subset the original Status
indicator to the value 1 (ie fully observed) or its survHE
rendition.
Now all the models are basically in alignment (as they should!).
I think, more interestingly, the point is that flexsurv
can handle the Status
indicator when you subset it to be =1 while still considering the fully observed as well as the censored cases. In that respect, flexsurv
is more general (and probably better) than survHE
. But I think the specificity of survHE
should be maintained --- even if at the expense of asking the user to format the data in the "correct" way...
Hi Gianluca
When I run:
m1<-fit.models(formula, data=regdata, distr="weibull")
plot(m1)
everything looks fine. However,m2<-fit.models(formula, data=regdata, distr="weibull", method="hmc")
plot(m2)
The model survival curves don't coincide at all with the KM as below. Is this a glitch or do I need to do something different?