Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

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

Open
helmingstay opened this issue Aug 9, 2024 · 3 comments
Open

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

helmingstay opened this issue Aug 9, 2024 · 3 comments

Comments

@helmingstay
Copy link

helmingstay commented Aug 9, 2024

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
Copy link
Owner

@mikesweeting Would you mind looking at this please?

@mikesweeting
Copy link
Contributor

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
Copy link
Owner

Thanks @mikesweeting - merged #196

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

3 participants