hesim-dev / hesim

Health economic simulation modeling and decision analysis
https://hesim-dev.github.io/hesim/
62 stars 15 forks source link

Allow for survival models with flexible baseline hazards and frailties #49

Closed canuckafar closed 3 years ago

canuckafar commented 3 years ago

The flexsurv package used for the survival models in hesim has many nice features, including the flexsurvspline function that allow for flexible fit to the hazard. However, it does not allow for frailties, which are important when modeling multilevel survival data, such as data from multiple different clinical trials. The stan_surv function of the rstanarm package can fit models with both a flexible baseline hazard and frailties. Could hesim incorporate this functionality?

dincerti commented 3 years ago

Hi, hesim could incorporate this functionality. The difficulty would depend on the baseline hazard. If splines (e.g., M-splines, B-splines) were required, then it might take a little more work. Also, it doesn't look like stan_surv() is part of the master branch in rstanarm yet?

If you wanted to use a parametric baseline hazard (e.g., Weibull), then you might be able to input the parameter values manually into a params_surv object using the current package functionality. For instance, if the frailty term was capturing trial heterogeneity, then (I think) the frailty term could be absorbed into the intercept.

canuckafar commented 3 years ago

First, thank you for your response.

You are correct, stan_surv is still developmental. But it works very nicely and can fit models very hard to fit otherwise like a model with M-spline or B-spline baseline hazards (which are very useful in cancer models where the baseline hazard can be very time dependent) and a hierarchical structure (i.e. frailties). What about the following work-around?

Assume a partitioned survival model, similar to the one in your package vignette for example, except with the following differences.

1) I fit a Bayesian survival model using some other R package, for each of the endpoints and each of the policies/strategies to be tested. Let's assume it is a cancer-related decision model, that there are two strategies to compare (chemotherapy vs immunotherapy) and two models (overall survival and progression-free survival) for each strategy (i.e. a simple stable-progressed-dead 3-state model) from which I will partition time. This means 4 survival models: chemo-OS, chemo-PFS, immuno-OS, and immuno-PFS.

2) I take the fitted survival models make a large number of draws from the posterior for each of the 4 survival fits, say 10,000 individual survival curves each, and save these 10,000 x 4 = 40,000 predicted survival curves in a data structure in R (I like grouped/nested tibbles for example).

3) I model the costs and utilities of the stable and progressed health states (maybe I get fancy and add a death cost to the dead state) using typical methods with average costs and utilities following some informed and reasonable distribution (maybe gamma and normal, respectively).

4) Is there a way to rig hesim so that I can get my simulated survival models into a Psm_Curves object so that I can merge it with the State_values objects containing my costs and utiliities?

Brant

dincerti commented 3 years ago

Hi Brant,

Actually that can work with a really small tweak which I've pushed to the stateprobs branch and will add to the development branch once I've had a time to properly test it out.

Per hesim 0.5.0, you can store your predicted survival curves in a survival object. I then created a sim_stateprobs.survival() function in the new branch which will compute state probabilities from the survival curves (like Psm$sim_stateprobs() currently does). You can then simply add the resulting stateprobs object to the Psm class object. Here is an example:

library("data.table")
hesim_dat <- hesim_data(
  strategies = data.table(strategy_id = 1:2L, 
                          strategy2 = c(0, 1)),
  patients = data.table(patient_id = 1L, cons = 1),
  states = data.table(state_id = 1:2L,
                      state_name = c("Stable", "Progression"))
)

# Generate survival curves, you can do this however you want
params1 <- params_surv(
  coefs = list(rate = as.matrix(data.frame(
    cons = c(log(.2), log(.25)),
    strategy2 = c(log(.8), log(.9))
  ))),
  dist = "exp"
)
params2 <- params_surv(
  coefs = list(rate = as.matrix(data.frame(
    cons = c(log(.1), log(.15)),
    strategy2 = c(log(.8), log(.9))
  ))),
  dist = "exp"
)
params <- params_surv_list(params1, params2)
psm_curves <- create_PsmCurves(params, input_data = expand(hesim_dat))
surv <- psm_curves$survival(t = 0:20)

# Dummy utility model for example
utility_tbl <- stateval_tbl(
  data.table(state_id = 1:2,
             est = c(1, 1)
  ),
  dist = "fixed"
)
utilitymod <- create_StateVals(utility_tbl, 
                               hesim_data = hesim_dat,
                               n = 2)

# Compute state probabilities
stprobs <- sim_stateprobs(surv)

# Use Psm class to compute QALYs given state probabilities computed above
psm <- Psm$new(utility_model = utilitymod)
psm$stateprobs_ <- stprobs
psm$sim_qalys()
psm$qalys_
canuckafar commented 3 years ago

That sounds promising. I look forward to trying it out.

dincerti commented 3 years ago

@bainman the development version now supports this. You can construct a survival object with survival() and simulate state probabilities with sim_stateprobs.survival(). The sim_stateprobs() documentation provides an example of how to do this in the context of a partitioned survival analysis.

canuckafar commented 3 years ago

Terrific. I will check it out this week.

On Sat, Feb 27, 2021 at 6:28 PM Devin Incerti notifications@github.com wrote:

@bainman https://github.com/bainman the development version now supports this. You can construct a survival object with survival() https://hesim-dev.github.io/hesim/dev/reference/survival.html and simulate state probabilities with sim_stateprobs.survival() https://hesim-dev.github.io/hesim/dev/reference/sim_stateprobs.survival.html. The sim_stateprobs() documentation provides an example of how to do this in the context of a partitioned survival analysis.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/hesim-dev/hesim/issues/49#issuecomment-787203570, or unsubscribe https://github.com/notifications/unsubscribe-auth/AS3DZTKPCWF5B7MECSHIQKDTBF5Y5ANCNFSM4XTOGBGQ .