aphp / heemod

Markov Models for Health Economic Evaluations
https://aphp.github.io/heemod/
GNU General Public License v3.0
11 stars 1 forks source link

Error with plot_values() function in heemod package #1

Open Avila-Diego opened 1 year ago

Avila-Diego commented 1 year ago

I am experiencing an error while using the plot_values() function from the heemod package in R. Specifically, when I attempt to use the panels = "by_value" argument, I receive the following error message:

> plot(result$model_runs, type = "values", panels = "by_state")
Error in plot.run_model(result$model_runs, type = "values", panels = "by_state") : 
  'panels' arguement not compatible.

Is there any possibility to tabulate or generate a graph that shows the cost and QALY by stage and by cycle of the entire simulation?

I am using heemod version 0.14.4 and R version 4.2.2.

KZARCA commented 1 year ago

Hi Diego,could you please send me a reproducible example? It will help me help you :)

Avila-Diego commented 1 year ago

Hello, thank you.

Let's suppose the following code:

library(heemod)
param <- define_parameters(
  age_init = 60,
  sex = 0,
  # age increases with cycles
  age = age_init + markov_cycle,

  # operative mortality rates
  omrPTHR = .02,
  omrRTHR = .02,

  # re-revision mortality rate
  rrr = .04,

  # parameters for calculating primary revision rate
  cons = -5.49094,
  ageC = -.0367,
  maleC = .768536,
  lambda = exp(cons + ageC * age_init + maleC * sex),
  gamma = 1.45367786,

  rrNP1 = .260677,

  # revision probability of primary procedure
  standardRR = 1 - exp(lambda * ((markov_cycle - 1) ^ gamma -
                                   markov_cycle ^ gamma)),
  np1RR = 1 - exp(lambda * rrNP1 * ((markov_cycle - 1) ^ gamma - 
                                      markov_cycle ^ gamma)),

  # age-related mortality rate
  sex_cat = ifelse(sex == 0, "FMLE", "MLE"),
  mr = get_who_mr(age, sex_cat,
                  country = "GBR", local = TRUE),

  # state values
  u_SuccessP = .85,
  u_RevisionTHR = .30,
  u_SuccessR = .75,
  c_RevisionTHR = 5294
)

mat_standard <- define_transition(
  state_names = c(
    "PrimaryTHR",
    "SuccessP",
    "RevisionTHR",
    "SuccessR",
    "Death"
  ),
  0, C, 0,          0, omrPTHR,
  0, C, standardRR, 0, mr,
  0, 0, 0,          C, omrRTHR+mr,
  0, 0, rrr,        C, mr,
  0, 0, 0,          0, 1
)

mat_np1 <- define_transition(
  state_names = c(
    "PrimaryTHR",
    "SuccessP",
    "RevisionTHR",
    "SuccessR",
    "Death"
  ),
  0, C, 0,          0, omrPTHR,
  0, C, np1RR,      0, mr,
  0, 0, 0,          C, omrRTHR+mr,
  0, 0, rrr,        C, mr,
  0, 0, 0,          0, 1
)

mod_standard <- define_strategy(
  transition = mat_standard,
  PrimaryTHR = define_state(
    utility = 0,
    cost = 394
  ),
  SuccessP = define_state(
    utility = discount(u_SuccessP, .015),
    cost = 0
  ),
  RevisionTHR = define_state(
    utility = discount(u_RevisionTHR, .015),
    cost = discount(c_RevisionTHR, .06)
  ),
  SuccessR = define_state(
    utility = discount(u_SuccessR, .015),
    cost = 0
  ),
  Death = define_state(
    utility = 0,
    cost = 0
  )
)

mod_np1 <- define_strategy(
  transition = mat_np1,
  PrimaryTHR = define_state(
    utility = 0,
    cost = 579
  ),
  SuccessP = define_state(
    utility = discount(u_SuccessP, .015),
    cost = 0
  ),
  RevisionTHR = define_state(
    utility = discount(u_RevisionTHR, .015),
    cost = discount(c_RevisionTHR, .06)
  ),
  SuccessR = define_state(
    utility = discount(u_SuccessR, .015),
    cost = 0
  ),
  Death = define_state(
    utility = 0,
    cost = 0
  )
)

res_mod <- run_model(
  standard = mod_standard,
  np1 = mod_np1,
  parameters = param,
  cycles = 60,
  cost = cost,
  effect = utility,
  method = "beginning"
)

The following graphics work well

plot(res_mod, panels = "by_state")
plot(res_mod, type = "values")

But when I want to obtain the costs and profits for each of the stadiums, I get an error

plot(res_mod$model_runs, type = "values", panels = "by_state")

I am sure it is not implemented, if you could guide me to get it, maybe I could develop it from the data already stored in memory.

KZARCA commented 1 year ago

I get it. I'm quite busy right now, and cannot code this function. If you would like to contribute, please check the function plot.run_model within R/strategy_print.R