chjackson / flexsurv

The flexsurv R package for flexible parametric survival and multi-state modelling
http://chjackson.github.io/flexsurv/
55 stars 27 forks source link

standsurv fails with bootstrap=T and cl=0.99 #195

Open helmingstay opened 3 months ago

helmingstay commented 3 months ago

I ran into the following error trying to compute 99% CI for model contrasts with standsurv():

Calculating bootstrap standard errors / confidence intervals
Error in `rename()`:
! Can't rename columns that don't exist.
✖ Column `2.5%` doesn't exist.

I built a minimal(ish) example from examples in the standsurv vignette [1]:

data(bc)
set.seed(236236)
## Age at diagnosis is correlated with survival time. A longer survival time 
## gives a younger mean age
bc$age <- rnorm(dim(bc)[1], mean = 65 - scale(bc$recyrs, scale=F), sd = 5)
## Create age at diagnosis in days - used later for matching to expected rates
bc$agedays <- floor(bc$age * 365.25)
## Create some random diagnosis dates between 01/01/1984 and 31/12/1989
bc$diag <- as.Date(floor(runif(dim(bc)[1], as.Date("01/01/1984", "%d/%m/%Y"), 
                               as.Date("31/12/1989", "%d/%m/%Y"))), 
                   origin="1970-01-01")
## Create sex (assume all are female)
bc$sex <- factor("female")
## 2-level prognostic variable
bc$group2 <- ifelse(bc$group=="Good", "Good", "Medium/Poor")

mod <- flexsurvreg(
    Surv(recyrs, censrec)~group2, 
    #anc = list(shape = ~ group2), dist="weibullPH"
    dist="weibull", data=bc
)

## example works either if commenting out line `B=...` or using `cl=0.95`
mod.contr <- standsurv(mod, 
        at=list(list(group2='Good'), list(group2='Medium/Poor')),
        B=1e4, boot=T,
        t=20, cl=0.99, ci=T 
)

Related: My original question was focused on family-wise error rates for model contrasts, similar to adjust="tukey" in emmeans() [2]. It seemed like a very simple adjustment would be to use wider CI (akin to a Bonferroni correction). I'd be interested to hear anyone's advice on using standsurv() for formal hypothesis testing (e.g. p-value of hazard ratio relative to 1).

[1] https://cran.r-project.org/web/packages/flexsurv/vignettes/standsurv.html [2] https://cran.r-project.org/web/packages/emmeans/vignettes/comparisons.html#pairwise

chjackson commented 3 months ago

@mikesweeting Would you mind looking at this please?

mikesweeting commented 3 months ago

Thank you @helmingstay for reporting this bug. I've now made the fix and created a pull request for @chjackson to incorporate. Thanks for spotting this!

I've not considered using standsurv for formal hypothesis testing so cannot help there.

chjackson commented 3 months ago

Thanks @mikesweeting - merged https://github.com/chjackson/flexsurv/pull/196