Error: can't find formula when estimating frailty model using stan_surv

[jburos]

Summary:

I am hitting an error when working with a shared-frailty model fit using stan_surv:

Error: Can't find formula stored in model frame or call

Description:

The formula isn't correctly recovered from the fit object, causing a number of different posterior operations to fail.

It appears to be caused in this line of formula.stanreg:

formula.stanreg <- function(x, ..., m = NULL) {
  if (is.mer(x) && !isTRUE(x$stan_function == "stan_gamm4")) 
      return(formula_mer(x, ...))   # <----  HERE
  if (is.surv(x))
    return(formula_surv(x, ...))
  x$formula
}

When the survival model contains random effects, the returned fit object is both is.mer & is.surv (both are true). The function should either reverse the order of the checks or exclude the stan_surv case as is done for stan_gamm4.

> # from example below:
> rstanarm:::is.mer(fit1)
[1] TRUE

Reproducible Steps:

library(survival)
library(tidyverse)
library(rstanarm) # from feature/survival branch
#> Loading required package: Rcpp
#> This is rstanarm version 2.21.2
#> - See https://mc-stan.org/rstanarm/articles/priors for changes to default priors!
#> - Default priors may change, so it's safest to specify priors, even if equivalent to the defaults.
#> - For execution on a local, multicore CPU with excess RAM we recommend calling
#>   options(mc.cores = parallel::detectCores())

data(leukemia, package = "survival")
#> Warning in data(leukemia, package = "survival"): data set 'leukemia' not found
leukemia <- tibble(leukemia) %>%
  mutate(id = seq_len(n())) %>%
  select(id, everything())

# fit model (works)
fit1 <- stan_surv(Surv(time = time, event = status) ~ x + (1 | id), data = leukemia)
#> 
#> SAMPLING FOR MODEL 'surv' NOW (CHAIN 1).
#> Chain 1: 
#> Chain 1: Gradient evaluation took 0.000117 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 1.17 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1: 
#> Chain 1: 
#> Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
#> Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
#> Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
#> Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
#> Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
#> Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
#> Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
#> Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
#> Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
#> Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
#> Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
#> Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
#> Chain 1: 
#> Chain 1:  Elapsed Time: 0.440674 seconds (Warm-up)
#> Chain 1:                0.391088 seconds (Sampling)
#> Chain 1:                0.831762 seconds (Total)
#> Chain 1: 
#> 
#> SAMPLING FOR MODEL 'surv' NOW (CHAIN 2).
#> Chain 2: 
#> Chain 2: Gradient evaluation took 2.4e-05 seconds
#> Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.24 seconds.
#> Chain 2: Adjust your expectations accordingly!
#> Chain 2: 
#> Chain 2: 
#> Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)
#> Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
#> Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
#> Chain 2: Iteration:  600 / 2000 [ 30%]  (Warmup)
#> Chain 2: Iteration:  800 / 2000 [ 40%]  (Warmup)
#> Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
#> Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
#> Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
#> Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
#> Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
#> Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
#> Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
#> Chain 2: 
#> Chain 2:  Elapsed Time: 0.390542 seconds (Warm-up)
#> Chain 2:                0.611345 seconds (Sampling)
#> Chain 2:                1.00189 seconds (Total)
#> Chain 2: 
#> 
#> SAMPLING FOR MODEL 'surv' NOW (CHAIN 3).
#> Chain 3: 
#> Chain 3: Gradient evaluation took 2.5e-05 seconds
#> Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.25 seconds.
#> Chain 3: Adjust your expectations accordingly!
#> Chain 3: 
#> Chain 3: 
#> Chain 3: Iteration:    1 / 2000 [  0%]  (Warmup)
#> Chain 3: Iteration:  200 / 2000 [ 10%]  (Warmup)
#> Chain 3: Iteration:  400 / 2000 [ 20%]  (Warmup)
#> Chain 3: Iteration:  600 / 2000 [ 30%]  (Warmup)
#> Chain 3: Iteration:  800 / 2000 [ 40%]  (Warmup)
#> Chain 3: Iteration: 1000 / 2000 [ 50%]  (Warmup)
#> Chain 3: Iteration: 1001 / 2000 [ 50%]  (Sampling)
#> Chain 3: Iteration: 1200 / 2000 [ 60%]  (Sampling)
#> Chain 3: Iteration: 1400 / 2000 [ 70%]  (Sampling)
#> Chain 3: Iteration: 1600 / 2000 [ 80%]  (Sampling)
#> Chain 3: Iteration: 1800 / 2000 [ 90%]  (Sampling)
#> Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
#> Chain 3: 
#> Chain 3:  Elapsed Time: 0.349901 seconds (Warm-up)
#> Chain 3:                0.409964 seconds (Sampling)
#> Chain 3:                0.759865 seconds (Total)
#> Chain 3: 
#> 
#> SAMPLING FOR MODEL 'surv' NOW (CHAIN 4).
#> Chain 4: 
#> Chain 4: Gradient evaluation took 2.5e-05 seconds
#> Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.25 seconds.
#> Chain 4: Adjust your expectations accordingly!
#> Chain 4: 
#> Chain 4: 
#> Chain 4: Iteration:    1 / 2000 [  0%]  (Warmup)
#> Chain 4: Iteration:  200 / 2000 [ 10%]  (Warmup)
#> Chain 4: Iteration:  400 / 2000 [ 20%]  (Warmup)
#> Chain 4: Iteration:  600 / 2000 [ 30%]  (Warmup)
#> Chain 4: Iteration:  800 / 2000 [ 40%]  (Warmup)
#> Chain 4: Iteration: 1000 / 2000 [ 50%]  (Warmup)
#> Chain 4: Iteration: 1001 / 2000 [ 50%]  (Sampling)
#> Chain 4: Iteration: 1200 / 2000 [ 60%]  (Sampling)
#> Chain 4: Iteration: 1400 / 2000 [ 70%]  (Sampling)
#> Chain 4: Iteration: 1600 / 2000 [ 80%]  (Sampling)
#> Chain 4: Iteration: 1800 / 2000 [ 90%]  (Sampling)
#> Chain 4: Iteration: 2000 / 2000 [100%]  (Sampling)
#> Chain 4: 
#> Chain 4:  Elapsed Time: 0.393629 seconds (Warm-up)
#> Chain 4:                0.337407 seconds (Sampling)
#> Chain 4:                0.731036 seconds (Total)
#> Chain 4:

# print fitted model (error)
fit1
#> stan_surv
#>  baseline hazard: M-splines on hazard scale
#> Error: Can't find formula stored in model frame or call.

# use posterior_survfit (error)
ppsurv <- posterior_survfit(fit1)
#> Error in terms.formula(formula, data = data): argument is not a valid model

# source of error:
formula(fit1)
#> Error: Can't find formula stored in model frame or call.

^{Created on 2021-06-17 by the reprex package (v2.0.0)}

RStanARM Version:

2.21.2, installed from github (stan-dev/rstanarm@dab004f)

R Version:

 setting  value                       
 version  R version 4.0.3 (2020-10-10)
 os       macOS Catalina 10.15.7      
 system   x86_64, darwin17.0          
 ui       RStudio                     
 language (EN)                        
 collate  en_US.UTF-8                 
 ctype    en_US.UTF-8                 
 tz       America/New_York            
 date     2021-06-17         

Operating System:

macOS Catalina 10.15.7

[jgabry]

Hey Jacki, thanks for the reproducible example and tracking down the source of the error!

The function should either reverse the order of the checks or exclude the stan_surv case as is done for stan_gamm4.

Yeah I agree that either of those options should work. @sambrilleman any preference or alternative solution?

[jburos]

Hi Jonah and Sam,

It looks like this pattern is repeated across a number of methods for stanreg objects (terms, for one). I'm stepping through a few of them to try to identify the set that needs to be updated.

Jacki

On Thu, Jun 17, 2021, 5:39 PM Jonah Gabry @.***> wrote:

Hey Jacki, thanks for the reproducible example and tracking down the source of the error!

The function should either reverse the order of the checks or exclude the stan_surv case as is done for stan_gamm4.

Yeah I think either of those options should work. @sambrilleman https://github.com/sambrilleman any preference or alternative solution?

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/stan-dev/rstanarm/issues/535#issuecomment-863581910, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAHBOPMFL5AZTJKWI4XGJMTTTJTSRANCNFSM464KXIEA .

[sambrilleman]

Hi Jacki

Awesome, thanks for picking up on the details of this! Makes it so much easier to fix, when you've essentially solved it! 😀

Is this a relatively recent install of stan_surv?

I think it has been introduced in this PR from a month ago: https://github.com/stan-dev/rstanarm/pull/526/files

That PR included a change to identify frailty models using is.mer() so that methods like ranef() would work. But it obviously had other consequences we didn't anticipate.

I'm happy to just switch the order of the statements if that works...?

[jburos]

Thanks @sambrilleman for the quick response! There are two places where the order needs to change, noted in this draft PR #536.

I will test how this impacts methods like ranef().

[jburos]

It looks like ranef() is failing on my test branch. I haven't looked into whether it worked prior to my changes.

[sambrilleman]

@jgabry Is it possible for us to turn on the automated CI tests for PRs to the feature/survival branch too (rather than just for PRs to master)? It might make it easier for us to catch some errors before deploying the stan_surv binary (I'm pretty sure it would have caught these ones...)

[sambrilleman]

I tried to take a closer look at this on the weekend, but I can't get a development version of rstanarm installed (and haven't been able to for a while). I didn't want to pollute this thread with those details, so I've opened a separate issue here: #537.

In the meantime, if @jburos can find a fix that works, then we can roll that out. Or if that is too messy then we can just roll back my earlier change to is.mer() in PR #526. Having posterior_survfit() etc work is more important than having ranef() work for now.

@jburos -- if the only thing appears to be that ranef() is broken, then that might be ok. I actually don't think I tested ranef() with the changes in PR #526, so maybe that fix never worked at all. 😦

[jburos]

Thanks @sambrilleman. Bummer, that stinks you can't get the branch to install. What a pain :/

I've been using my test branch (referenced in PR #536) without issues, other than the ranef function not working. I did step through that ranef function and encountered some other issues making me think it may not have worked previously. I can open a separate issue or add this detail to #440.

I'll run the tests locally for PR #536 & report the output there. That may help you decide whether to wait, merge this PR, or rollback the earlier change to is.mer. I don't have any preference there; just putting the PR up in case it's helpful.

[sambrilleman]

@jburos I think maybe we just merge your PR if it seems like everything is working except ranef. I think we are probably going to want/need is.mer() to work for frailty models at some point, so better to roll out your PR now which gets us some of the way there. We can figure out the fix for ranef in a later PR I guess, since it wasn't working anyway...!

[sambrilleman]

Also, is this bug in the binary hosted at https://mc-stan.org/r-packages/, or did you encounter this bug after installing directly from Github?

[jburos]

Thanks Sam. I installed directly from github, on the latest branch. I can check out the binary there; that looks helpful.

I can create an issue to restore/fix ranef behaviour.

On Tue, Jun 22, 2021 at 6:05 AM Sam Brilleman @.***> wrote:

Also, is this bug in the binary hosted at https://mc-stan.org/r-packages/, or did you encounter this bug after installing directly from Github?

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/stan-dev/rstanarm/issues/535#issuecomment-865849230, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAHBOPJGXXMAEF2ESYFFSSLTUBOAPANCNFSM464KXIEA .

[padpadpadpad]

Hi guys. Any progress on this?

I get these errors when running: remotes::install_github('stan-dev/rstanarm@feature/survival')

I tried to install the binary from https://mc-stan.org/r-packages/ using install.packages("rstanarm", repos = c("https://mc-stan.org/r-packages/", getOption("repos"))) then saying no when i did not want to compile the newer version from source. However this resulted in stan_surv not being present in the package.

Is there any legacy package that I could install to be able to run the models.

Cheers Dan

[sambrilleman]

Hi @padpadpadpad - yeah I think you are right - it seems that the version now hosted/installed via:

install.packages("rstanarm", repos = c("https://mc-stan.org/r-packages/", getOption("repos")))

is no longer the binary with stan_surv.

Any chance you know when that was updated/changed @jgabry or @rok-cesnovar?

Moving forward / as an alternative, do you think we are good to merge your PR @jburos? Based on your feedback from the tests, I'm happy to if you are. After that, at least @padpadpadpad could install using remotes::install_github('stan-dev/rstanarm@feature/survival') (at least until an updated binary with stan_surv is uploaded again).

[rok-cesnovar]

Any chance you know when that was updated/changed @jgabry or @rok-cesnovar?

Nothing has been changed, but it might pick up the one from CRAN as the versions there may have been bumped?

Try:

remove.packages("rstanarm")
install.packages("rstanarm", repos = c("https://mc-stan.org/r-packages/"))

The difference is this omits getOption("repos"), so it wont search CRAN.

[jburos]

Hi Sam,

Yes I think it's good to merge. I've been using it and haven't seen issues; also the problem with the frailty models is solved in that branch.

Going forward though I agree it would be great up set up some automated testing of this branch.

Jacki

On Sat, Jul 24, 2021, 8:19 PM Sam Brilleman @.***> wrote:

Hi @padpadpadpad https://github.com/padpadpadpad - yeah I think you are right - it seems that the version now hosted/installed via:

install.packages("rstanarm", repos = c("https://mc-stan.org/r-packages/", getOption("repos")))

is no longer the binary with stan_surv.

Any chance you know when that was updated/changed @jgabry https://github.com/jgabry or @rok-cesnovar https://github.com/rok-cesnovar?

Moving forward / as an alternative, do you think we are good to merge your PR https://github.com/stan-dev/rstanarm/pull/536 @jburos https://github.com/jburos? Based on your feedback from the tests, I'm happy to if you are. After that, at least @padpadpadpad https://github.com/padpadpadpad could install using @.***/survival') (at least until an updated binary with stan_surv is uploaded again).

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/stan-dev/rstanarm/issues/535#issuecomment-886127620, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAHBOPMXEK4SF3G2WH5QMTLTZNKC7ANCNFSM464KXIEA .

[jburos]

It's been almost a year now and I wanted to circle back to this. Is there any appetite for us to:

1. Merge the pending fixes to the frailty model estimation
2. Set up or turn on automatic testing for this survival branch?

I'm willing to donate some time to help on higher-priority fixes to this branch, as I think supporting Bayesian survival models is an important goal. I'm curious if there is interest from rstanarm more generally. Possibly automate the updating of the binaries (& documenting them better) as well?

cc @jgabry @rok-cesnovar @sambrilleman

[sambrilleman]

Hi @jburos - thanks for nudging this! Sorry I didn't just merge it earlier. I've just merged PR #536, which I guess would close this issue.

But in relation to the wider questions / comments:

Set up or turn on automatic testing for this survival branch? ... Possibly automate the updating of the binaries (& documenting them better) as well?

Yeah I agree - these would be great additions. I think @jgabry set up a CI pipeline for automated testing on the master branch, so would be great if we can extend it to the survival branch, and include publishing the binaries in that CI pipeline.

More broadly, I'd be keen to know what the view is amongst the Stan development and rstanarm team for stan_surv longer term. It is a bit sad that it has never made it onto CRAN for an official release.

@jgabry @bgoodri - are CRAN actually still blocking us from ever adding more models to rstanarm? Or is it more about a lack of appetite to take on the risk of adding new models and any bugs/maintenance that might introduce?

Realistically, I won't have much time to contribute to coding on this project moving forward. But I'm happy to be involved in some sort of advisory capacity, as it seems that @jburos and @ermeel86 (and perhaps others) also are. Maybe, through the Stan community, we can find someone (e.g. a Phd or postdoc) with more time on their hands and an interest/expertise in Bayesian survival analysis and R who could deal with maintenance or feature additions, and which might give us a way to reliably support adding these models to CRAN. It seems like a waste to just leave them to rust (although I guess they have been getting some use in their current hosted binary form, although I imagine it is a lot less exposure than if they were in the CRAN release proper).

But yeah if the CRAN option isn't realistic, then as @jburos mentioned, at a minimum it would be amazing if we could find a way to have a CI pipeline for the feature/survival branch that tests and publishes the Stan hosted binary. Perhaps we can try to find someone with the time and interest to work on that. If none of us, then maybe through the forums?

[bgoodri]

I think we could merge it now that we are utilizing LTO on all versions on Windows. I haven't tried it on 4.2.x but things somehow compile on all other versions of Windows. It is just that the LTO that is supported by CRAN does not cover C++, so we had to do some weird things in order to make it work. But CRAN has not objected yet. Also, 32bit is no longer a concern.

On Fri, Apr 29, 2022 at 1:51 AM Sam Brilleman @.***> wrote:

Hi @jburos https://github.com/jburos - thanks for nudging this! Sorry I didn't just merge it earlier. I've just merged PR #536 https://github.com/stan-dev/rstanarm/pull/536, which I guess would close this issue.

But in relation to the wider questions / comments:

Set up or turn on automatic testing for this survival branch? ... Possibly automate the updating of the binaries (& documenting them better) as well?

Yeah I agree - these would be great additions. I think @jgabry https://github.com/jgabry set up a CI pipeline for automated testing on the master branch, so would be great if we can extend it to the survival branch, and include publishing the binaries in that CI pipeline.

More broadly, I'd be keen to know what the view is amongst the Stan development and rstanarm team for stan_surv longer term. It is a bit sad that it has never made it onto CRAN for an official release.

@jgabry https://github.com/jgabry @bgoodri https://github.com/bgoodri

• are CRAN actually still blocking us from ever adding more models to rstanarm? Or is it more about a lack of appetite to take on the risk of adding new models and any bugs/maintenance that might introduce?

Realistically, I won't have much time to contribute to coding on this project moving forward. But I'm happy to be involved in some sort of advisory capacity, as it seems that @jburos https://github.com/jburos and @ermeel86 https://github.com/ermeel86 (and perhaps others) also are. Maybe, through the Stan community, we can find someone (e.g. a Phd or postdoc) with more time on their hands and an interest/expertise in Bayesian survival analysis and R who could deal with maintenance or feature additions, and which might give us a way to reliably support adding these models to CRAN. It seems like a waste to just leave them to rust (although I guess they have been getting some use in their current hosted binary form, although I imagine it is a lot less exposure than if they were in the CRAN release proper).

But yeah if the CRAN option isn't realistic, then as @jburos https://github.com/jburos mentioned, at a minimum it would be amazing if we could find a way to have a CI pipeline for the feature/survival branch that tests and publishes the Stan hosted binary. Perhaps we can try to find someone with the time and interest to work on that. If none of us, then maybe through the forums?

— Reply to this email directly, view it on GitHub https://github.com/stan-dev/rstanarm/issues/535#issuecomment-1112899327, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAZ2XKVOFH3JRQH6ZDNYZK3VHN2EZANCNFSM464KXIEA . You are receiving this because you were mentioned.Message ID: @.***>

[jburos]

This would be amazing. I'm happy to help if there are tests or other refactoring needed to get this merged into the master branch.

Maybe we should post on the stan discourse to see if there are others who can help both to get this into CRAN & with ongoing maintenance. There are a couple issues on github that we might want to address.

On Fri, Apr 29, 2022 at 9:37 AM bgoodri @.***> wrote:

I think we could merge it now that we are utilizing LTO on all versions on Windows. I haven't tried it on 4.2.x but things somehow compile on all other versions of Windows. It is just that the LTO that is supported by CRAN does not cover C++, so we had to do some weird things in order to make it work. But CRAN has not objected yet. Also, 32bit is no longer a concern.

On Fri, Apr 29, 2022 at 1:51 AM Sam Brilleman @.***> wrote:

Hi @jburos https://github.com/jburos - thanks for nudging this! Sorry I didn't just merge it earlier. I've just merged PR #536 https://github.com/stan-dev/rstanarm/pull/536, which I guess would close this issue.

But in relation to the wider questions / comments:

Set up or turn on automatic testing for this survival branch? ... Possibly automate the updating of the binaries (& documenting them better) as well?

Yeah I agree - these would be great additions. I think @jgabry https://github.com/jgabry set up a CI pipeline for automated testing on the master branch, so would be great if we can extend it to the survival branch, and include publishing the binaries in that CI pipeline.

More broadly, I'd be keen to know what the view is amongst the Stan development and rstanarm team for stan_surv longer term. It is a bit sad that it has never made it onto CRAN for an official release.

@jgabry https://github.com/jgabry @bgoodri <https://github.com/bgoodri

• are CRAN actually still blocking us from ever adding more models to rstanarm? Or is it more about a lack of appetite to take on the risk of adding new models and any bugs/maintenance that might introduce?

Realistically, I won't have much time to contribute to coding on this project moving forward. But I'm happy to be involved in some sort of advisory capacity, as it seems that @jburos https://github.com/jburos and @ermeel86 https://github.com/ermeel86 (and perhaps others) also are. Maybe, through the Stan community, we can find someone (e.g. a Phd or postdoc) with more time on their hands and an interest/expertise in Bayesian survival analysis and R who could deal with maintenance or feature additions, and which might give us a way to reliably support adding these models to CRAN. It seems like a waste to just leave them to rust (although I guess they have been getting some use in their current hosted binary form, although I imagine it is a lot less exposure than if they were in the CRAN release proper).

But yeah if the CRAN option isn't realistic, then as @jburos https://github.com/jburos mentioned, at a minimum it would be amazing if we could find a way to have a CI pipeline for the feature/survival branch that tests and publishes the Stan hosted binary. Perhaps we can try to find someone with the time and interest to work on that. If none of us, then maybe through the forums?

— Reply to this email directly, view it on GitHub <https://github.com/stan-dev/rstanarm/issues/535#issuecomment-1112899327 , or unsubscribe < https://github.com/notifications/unsubscribe-auth/AAZ2XKVOFH3JRQH6ZDNYZK3VHN2EZANCNFSM464KXIEA

. You are receiving this because you were mentioned.Message ID: @.***>

— Reply to this email directly, view it on GitHub https://github.com/stan-dev/rstanarm/issues/535#issuecomment-1113321142, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAHBOPNQNEU6KI5VAB3K4ADVHPQZ7ANCNFSM464KXIEA . You are receiving this because you were mentioned.Message ID: @.***>