search

r-lib / vctrs

Generic programming with typed R vectors
https://vctrs.r-lib.org
Other
282 stars 65 forks source link

Markov cost effectiveness model - survival curve is using a best fit line - I do not want it to. #1909

Closed mmaughan0 closed 5 months ago

mmaughan0 commented 5 months ago

I am building a time dependent markov cost effectiveness model comparing different ebola treatments over a 28 day period and have received the error above - I suspect the transition matrices is where the issue arises. I was able to reproduce this type of model from a paper, but am getting bugs when customizing for my own uses. (search "mat_base") if you want to go straight to the matrices.

Thanks!!!

library(heemod) library(flexsurv) library(rgho)

library(survminer)

par_mod <- define_parameters(

p_death_disease_base = compute_surv( fit_death_disease_base, time = state_time, km_limit = 28))

code below creates a variable in which mortality curves reside. This variable (p_death_disease_strat)

can then be used in transition probability matrices etc.

par_mod <- modify( par_mod,

p_death_disease_REGN_EB3 = compute_surv( fit_death_disease_REGNEB3, time = state_time, km_limit = 28))

par_mod <- modify( par_mod,

p_death_disease_mAb114 = compute_surv( fit_death_disease_mAb114, time = state_time, km_limit = 28))

pulling in induct to treat table

par_mod <- modify( par_mod,

p_tx_screen = compute_surv( fit_screen_to_treatment, time = state_time, km_limit = 28))

pulling in discharge tables

par_mod <- modify( par_mod,

p_discharge_base = compute_surv( fit_discharge_base, time = state_time, km_limit = 28))

par_mod <- modify( par_mod,

p_discharge_mAb114 = compute_surv( fit_discharge_mAb114, time = state_time, km_limit = 28))

par_mod <- modify( par_mod,

p_discharge_REGN_EB3 = compute_surv( fit_discharge_REGN_EB3, time = state_time, km_limit = 28))

This code is what imports WHO life tables

fit_death_disease_base <- flexsurv::flexsurvreg( survival::Surv(Time, Event) ~ 1, dist = "weibull", data = Base_Case_TTE_Table)

fit_death_disease_mAb114 <- flexsurv::flexsurvreg( survival::Surv(Time, Event) ~ 1, dist = "weibull", data = mAb114_TTE_Table)

fit_death_disease_REGNEB3 <- flexsurv::flexsurvreg( survival::Surv(Time, Event) ~ 1, dist = "weibull", data = REGN_EB3_TTE_Table)

fit_screen_to_treatment <- flexsurv::flexsurvreg( survival::Surv(Time, Event) ~ 1, dist = "weibull", data = Induct_to_Tx_TTE_table)

fit_discharge_base <- flexsurv::flexsurvreg( survival::Surv(Time, Event) ~ 1, dist = "weibull", data = Base_Case_Discharge_TTE_Table)

fit_discharge_mAb114 <- flexsurv::flexsurvreg( survival::Surv(Time, Event) ~ 1, dist = "weibull", data = mAb114_TTE_Discharge_Table)

fit_discharge_REGN_EB3 <- flexsurv::flexsurvreg( survival::Surv(Time, Event) ~ 1, dist = "weibull", data = mAb114_TTE_Discharge_Table)

plot(fit_discharge_REGN_EB3)

define matrices

putting in point estimates from Levine et. al. paper. All diagnostic rates can be derived from these 3 data points

sensitivity_algo = .93 specificity_algo = .23 prevalence = .42

positive_algo = sensitivity_algo prevalence + (1- specificity_algo) (1- prevalence)

p_pos_low = ((1-sensitivity_algo)prevalence)/(((1- sensitivity_algo) prevalence) + (specificity_algo (1- prevalence))) p_pos_hi = (sensitivity_algo prevalence)/((sensitivity_algo prevalence)+((1- prevalence)(1- specificity_algo)))

define base transition matrix

mat_base <- define_transition( state_names = c("induction_discharge", "triage_hi", "triage_lo", "conf_sick", "death"),

0, positive_algo, C, 0, 0, c, 0, 0, p_pos_hip_tx_screen, 0, c, 0, 0, p_pos_lowp_tx_screen, 0, p_discharge_base, 0, 0, c, p_death_disease_base, 0, 0, 0, 0, 1)

mat_mAb_114 <- define_transition( state_names = c("induction_discharge", "triage_hi", "triage_lo", "conf_sick", "death"),

0, positive_algo, C, 0, 0, 0, 0, C, p_pos_hip_tx_screen, 0, 0,0, C, p_pos_lowp_tx_screen, 0, p_discharge_mAb114, 0, 0, C, p_death_disease_mAb114, 0, 0, 0, 0, 1)

mat_REGN_EB3 <- define_transition( state_names = c("induction_discharge", "triage_hi", "triage_lo", "conf_sick", "death"),

0, positive_algo, C, 0, 0, 0, 0, C, p_pos_hip_tx_screen, 0, 0,0, C, p_pos_lowp_tx_screen, 0, p_discharge_REGN_EB3, 0, 0, C, p_death_disease_REGN_EB3, 0, 0, 0, 0, 1)

define state values

unknown represents induction - is a transient state that occurs between presentation and triage.

Patients do not spend a full cycle (day) there,so no qaly assigned

some marginal ppe usage assumed as well as clinician time that always has opportunity cost in this setting

state_induction_discharge <- define_state( cost_total = 0, qaly = 0 )

state_triage_hi <- define_state( cost_total = 150, qaly = 1/365 )

state_triage_low <- define_state( cost_total = 100, qaly = 1/365 )

Soc_cost = 200

state_conf_sick <- define_state( cost_total = Soc_cost, qaly = 1/365 )

before defining the sick state, we need to allow it to only charge for mAb on day 1

par_mod <- modify( par_mod,

cost_tx_start = 300, cost_tx_end = Soc_cost, n_days = 1, cost_tx_cycle_mAb114 = ifelse( state_time < n_days, cost_tx_start, cost_tx_end) )

state_sick_mAb114 <- define_state( cost_total = cost_tx_cycle_mAb114, qaly = 1/365 )

par_mod <- modify( par_mod,

cost_tx_start = 310, cost_tx_end = 200, n_days = 1.1, cost_tx_cycle_REGN_EB3 = ifelse( state_time < n_days, cost_tx_start, cost_tx_end) )

state_sick_REGN_Eb3 <- define_state( cost_total = cost_tx_cycle_REGN_EB3, qaly = 1/365 )

state_death <- define_state( cost_total = 50, qaly = -30 )

define strategies

strat_base <- define_strategy( transition = mat_base,

induction_discharge = state_induction_discharge, triage_hi = state_triage_hi, triage_lo = state_triage_low, conf_sick = state_conf_sick, death = state_death )

strat_mAb114 <- define_strategy( transition = mat_mAb_114,

induction_discharge = state_induction_discharge, triage_hi = state_triage_hi, triage_lo = state_triage_low, conf_sick = state_sick_mAb114, death = state_death )

strat_REGN_EB3 <- define_strategy( transition = mat_REGN_EB3,

induction_discharge = state_induction_discharge, triage_hi = state_triage_hi, triage_lo = state_triage_low, conf_sick = state_sick_REGN_Eb3, death = state_death )

res_mod <- run_model( parameters = par_mod,

mAb114 = strat_mAb114, REGN_EB3 = strat_REGN_EB3, base = strat_base,

cycles = 30,

cost = cost_total, effect = qaly,

method = "life-table" )

DavisVaughan commented 5 months ago

Could you please turn this into a self-contained reprex (short for minimal reproducible example)? It will help us help you if we can be sure we're all working with/looking at the same stuff.

If you've never heard of a reprex before, you might want to start by reading the tidyverse.org help page.

You can install reprex by running (you may already have it, though, if you have the tidyverse package installed):

install.packages("reprex")

Thanks

jennybc commented 5 months ago

Cmd/Ctrl + V?

mmaughan0 commented 5 months ago

Ok it works now

library(heemod)
library(flexsurv)
#> Loading required package: survival
library(rgho)

objects
#> function (name, pos = -1L, envir = as.environment(pos), all.names = FALSE, 
#>     pattern, sorted = TRUE) 
#> {
#>     if (!missing(name)) {
#>         pos <- tryCatch(name, error = function(e) e)
#>         if (inherits(pos, "error")) {
#>             name <- substitute(name)
#>             if (!is.character(name)) 
#>                 name <- deparse(name)
#>             warning(gettextf("%s converted to character string", 
#>                 sQuote(name)), domain = NA)
#>             pos <- name
#>         }
#>     }
#>     all.names <- .Internal(ls(envir, all.names, sorted))
#>     if (!missing(pattern)) {
#>         if ((ll <- length(grep("[", pattern, fixed = TRUE))) && 
#>             ll != length(grep("]", pattern, fixed = TRUE))) {
#>             if (pattern == "[") {
#>                 pattern <- "\\["
#>                 warning("replaced regular expression pattern '[' by  '\\\\['")
#>             }
#>             else if (length(grep("[^\\\\]\\[<-", pattern))) {
#>                 pattern <- sub("\\[<-", "\\\\\\[<-", pattern)
#>                 warning("replaced '[<-' by '\\\\[<-' in regular expression pattern")
#>             }
#>         }
#>         grep(pattern, all.names, value = TRUE)
#>     }
#>     else all.names
#> }
#> <bytecode: 0x000001da484a5dd0>
#> <environment: namespace:base>

##par_mod <- define_parameters(
  ##age_base = 30*365,##cycle in this model is 1 day
  ##age_cycle = model_time + age_base)

par_mod <- define_parameters(

  p_death_disease_base = compute_surv(
    fit_death_disease_base,
    time = state_time,
    km_limit = 28))

##code below creates a variable in which mortality curves reside. This variable (p_death_disease_strat) 
##can then be used in transition probability matrices etc.##

par_mod <- modify(
  par_mod,

  p_death_disease_REGN_EB3 = compute_surv(
    fit_death_disease_REGNEB3,
    time = state_time,
    km_limit = 28))

par_mod <- modify(
  par_mod,

  p_death_disease_mAb114 = compute_surv(
     fit_death_disease_mAb114,
     time = state_time,
     km_limit = 28))

##pulling in induct to treat table##

par_mod <- modify(
  par_mod,

  p_tx_screen = compute_surv(
    fit_screen_to_treatment,
    time = state_time,
    km_limit = 28))

###pulling in discharge tables###

par_mod <- modify(
  par_mod,

  p_discharge_base = compute_surv(
    fit_discharge_base,
    time = state_time,
    km_limit = 28))

par_mod <- modify(
  par_mod,

  p_discharge_mAb114 = compute_surv(
    fit_discharge_mAb114,
    time = state_time,
    km_limit = 28))

par_mod <- modify(
  par_mod,

  p_discharge_REGN_EB3 = compute_surv(
    fit_discharge_REGN_EB3,
    time = state_time,
    km_limit = 28))

##This code is what imports WHO life tables##

fit_death_disease_base <- flexsurv::flexsurvreg(
  survival::Surv(Time, Event) ~ 1,
  dist = "weibull",
  data = Base_Case_TTE_Table)
#> Error in eval(expr, envir, enclos): object 'Base_Case_TTE_Table' not found

fit_death_disease_mAb114 <- flexsurv::flexsurvreg(
  survival::Surv(Time, Event) ~ 1,
  dist = "weibull",
  data = mAb114_TTE_Table)
#> Error in eval(expr, envir, enclos): object 'mAb114_TTE_Table' not found

fit_death_disease_REGNEB3 <- flexsurv::flexsurvreg(
  survival::Surv(Time, Event) ~ 1,
  dist = "weibull",
  data = REGN_EB3_TTE_Table)
#> Error in eval(expr, envir, enclos): object 'REGN_EB3_TTE_Table' not found

fit_screen_to_treatment <- flexsurv::flexsurvreg(
  survival::Surv(Time, Event) ~ 1,
  dist = "weibull",
  data = Induct_to_Tx_TTE_table)
#> Error in eval(expr, envir, enclos): object 'Induct_to_Tx_TTE_table' not found

fit_discharge_base <- flexsurv::flexsurvreg(
  survival::Surv(Time, Event) ~ 1,
  dist = "weibull",
  data = Base_Case_Discharge_TTE_Table)
#> Error in eval(expr, envir, enclos): object 'Base_Case_Discharge_TTE_Table' not found

fit_discharge_mAb114 <- flexsurv::flexsurvreg(
  survival::Surv(Time, Event) ~ 1,
  dist = "weibull",
  data = mAb114_TTE_Discharge_Table)
#> Error in eval(expr, envir, enclos): object 'mAb114_TTE_Discharge_Table' not found

fit_discharge_REGN_EB3 <- flexsurv::flexsurvreg(
  survival::Surv(Time, Event) ~ 1,
  dist = "weibull",
  data = mAb114_TTE_Discharge_Table)
#> Error in eval(expr, envir, enclos): object 'mAb114_TTE_Discharge_Table' not found

###define matrices####

##putting in point estimates from Levine et. al. paper. All diagnostic rates can be derived from these 3 data points##

sensitivity_algo = .93
specificity_algo = .23
prevalence = .42

positive_algo = sensitivity_algo * prevalence + (1- specificity_algo) * (1- prevalence)

p_pos_low = ((1-sensitivity_algo)*prevalence)/(((1- sensitivity_algo)* prevalence) + (specificity_algo * (1- prevalence)))
p_pos_hi = (sensitivity_algo * prevalence)/((sensitivity_algo * prevalence)+((1- prevalence)*(1- specificity_algo)))

##define base transition matrix##

mat_base <- define_transition(
  state_names = c("induction_discharge", "triage_hi", "triage_lo", "conf_sick", "death"),

  0, positive_algo, C, 0, 0,
  c, 0, 0, p_pos_hi*p_tx_screen, 0, 
  c, 0, 0, p_pos_low*p_tx_screen, 0, 
  p_discharge_base, 0, 0, c, p_death_disease_base,
  0, 0, 0, 0, 1)

mat_mAb_114 <- define_transition(
  state_names = c("induction_discharge", "triage_hi", "triage_lo", "conf_sick", "death"),

  0, positive_algo, C, 0, 0,
  0, 0, C, p_pos_hi*p_tx_screen, 0,
  0,0, C, p_pos_low*p_tx_screen, 0,
  p_discharge_mAb114, 0, 0, C, p_death_disease_mAb114,
  0, 0, 0, 0, 1)

mat_REGN_EB3 <- define_transition(
  state_names = c("induction_discharge", "triage_hi", "triage_lo", "conf_sick", "death"),

  0, positive_algo, C, 0, 0,
  0, 0, C, p_pos_hi*p_tx_screen, 0,
  0,0, C, p_pos_low*p_tx_screen, 0,
  p_discharge_REGN_EB3, 0, 0, C, p_death_disease_REGN_EB3,
  0, 0, 0, 0, 1)

###define state values###

##unknown represents induction - is a transient state that occurs between presentation and triage. 
##Patients do not spend a full cycle (day) there,so no qaly assigned
##some marginal ppe usage assumed as well as clinician time that always has opportunity cost in this setting

state_induction_discharge <- define_state(
  cost_total = 0, 
  qaly = 0
  )

state_triage_hi <- define_state(
  cost_total = 150,
  qaly = 1/365
)

state_triage_low <- define_state(
  cost_total = 100,
  qaly = 1/365
)

Soc_cost = 200 

state_conf_sick <- define_state(
 cost_total = Soc_cost, 
  qaly = 1/365
)

###before defining the sick state, we need to allow it to only charge for mAb on day 1###

par_mod <- modify(
  par_mod,

  cost_tx_start = 300,
  cost_tx_end = Soc_cost,
  n_days = 1,
  cost_tx_cycle_mAb114 = ifelse(
    state_time < n_days,
    cost_tx_start,
    cost_tx_end)
  )

state_sick_mAb114 <- define_state(
  cost_total = cost_tx_cycle_mAb114,
  qaly = 1/365
)

par_mod <- modify(
  par_mod,

  cost_tx_start = 310,
  cost_tx_end = 200,
  n_days = 1.1,
  cost_tx_cycle_REGN_EB3 = ifelse(
    state_time < n_days,
    cost_tx_start,
    cost_tx_end)
)

state_sick_REGN_Eb3 <- define_state(
  cost_total = cost_tx_cycle_REGN_EB3,
  qaly = 1/365
)

state_death <- define_state(
  cost_total = 50,
  qaly = -30
)

##define strategies##

strat_base <- define_strategy(
  transition = mat_base,

  induction_discharge = state_induction_discharge,
  triage_hi = state_triage_hi,
  triage_lo = state_triage_low,
  conf_sick = state_conf_sick,
  death = state_death
)

strat_mAb114 <- define_strategy(
  transition = mat_mAb_114,

  induction_discharge = state_induction_discharge,
  triage_hi = state_triage_hi,
  triage_lo = state_triage_low,
  conf_sick = state_sick_mAb114,
  death = state_death
)

strat_REGN_EB3 <- define_strategy(
  transition = mat_REGN_EB3,

  induction_discharge = state_induction_discharge,
  triage_hi = state_triage_hi,
  triage_lo = state_triage_low,
  conf_sick = state_sick_REGN_Eb3,
  death = state_death
)

res_mod <- run_model(
  parameters = par_mod,

  mAb114 = strat_mAb114,
  REGN_EB3 = strat_REGN_EB3,
  base = strat_base,

  cycles = 30,

  cost = cost_total,
  effect = qaly,

  method = "life-table"
)
#> mAb114: detected use of 'state_time', expanding states: conf_sick, triage_hi, triage_lo.
#> Error: Error in parameter: p_death_disease_base: Error in dplyr::mutate(start_tibble, !!!x_tidy[seq_len(i)]) : 
#>   ℹ In argument: `p_death_disease_base =
#>   compute_surv(fit_death_disease_base, time = state_time, km_limit = 28)`.
#> Caused by error in `FUN()`:
#> ! object 'fit_death_disease_base' not found

Created on 2024-01-17 with reprex v2.1.0

mmaughan0 commented 5 months ago

also I am not getting the errors when I pull in the tables on my end FYI - my code runs fine except for the final error as mentioned above "Error in rep(res, nr) : attempt to replicate an object of type 'builtin'"

THanks I see how this is way more readable

mmaughan0 commented 5 months ago

Cmd/Ctrl + V?

hahahah I am not THAT much of a novice. I didnt see the preview pane and how it showed in the reprex format. I thought I had just copied and pasted the code like I had done in my original comment because I did not see the preview pane. Thanks!!!

DavisVaughan commented 5 months ago

Do you see how the reprex output errors with "object not found" errors? That means you haven't included enough to be able to run that code in a standalone way. You may need to manually create some "fake" data frames to be able to run it, or you could try and figure out the exact function that seems to be the issue and figure out what you are giving it as input.

Without a working reprex, it is fairly hard for us to figure out how to help!

mmaughan0 commented 5 months ago

Understood. I pulled in excel tables for my real example but I can mock up dummy tables in the code. will get another reprex and shoot it back.

mmaughan0 commented 5 months ago 
library(heemod)
library(flexsurv)
#> Loading required package: survival
library(rgho)

objects
#> function (name, pos = -1L, envir = as.environment(pos), all.names = FALSE, 
#>     pattern, sorted = TRUE) 
#> {
#>     if (!missing(name)) {
#>         pos <- tryCatch(name, error = function(e) e)
#>         if (inherits(pos, "error")) {
#>             name <- substitute(name)
#>             if (!is.character(name)) 
#>                 name <- deparse(name)
#>             warning(gettextf("%s converted to character string", 
#>                 sQuote(name)), domain = NA)
#>             pos <- name
#>         }
#>     }
#>     all.names <- .Internal(ls(envir, all.names, sorted))
#>     if (!missing(pattern)) {
#>         if ((ll <- length(grep("[", pattern, fixed = TRUE))) && 
#>             ll != length(grep("]", pattern, fixed = TRUE))) {
#>             if (pattern == "[") {
#>                 pattern <- "\\["
#>                 warning("replaced regular expression pattern '[' by  '\\\\['")
#>             }
#>             else if (length(grep("[^\\\\]\\[<-", pattern))) {
#>                 pattern <- sub("\\[<-", "\\\\\\[<-", pattern)
#>                 warning("replaced '[<-' by '\\\\[<-' in regular expression pattern")
#>             }
#>         }
#>         grep(pattern, all.names, value = TRUE)
#>     }
#>     else all.names
#> }
#> <bytecode: 0x000002356f1c75a0>
#> <environment: namespace:base>

##par_mod <- define_parameters(
##age_base = 30*365,##cycle in this model is 1 day
##age_cycle = model_time + age_base)

par_mod <- define_parameters(

  p_death_disease_base = compute_surv(
    fit_death_disease_base,
    time = state_time,
    km_limit = 14))

#establishing dummy table only for purpose of reprex. In reality each strategy would have a different table

tab_surv <- structure(list(time = c(0.4, 8.7, 7, 5.1, 9.2, 1, 0.5, 3.3, 1.8,
                                    3, 6.7, 3.7, 1.1, 5.9, 5.1, 10, 10, 10, 10, 10, 10, 10, 10, 10,
                                    10), status = c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
                                                    1L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L)), .Names = c("time",
                                                                                                                     "status"), row.names = c(NA, -25L), class = "data.frame")

##code below creates a variable in which mortality curves reside. This variable (p_death_disease_strat) 
##can then be used in transition probability matrices etc.##

par_mod <- modify(
  par_mod,

  p_death_disease_REGN_EB3 = compute_surv(
    fit_death_disease_REGNEB3,
    time = state_time,
    km_limit = 28))

par_mod <- modify(
  par_mod,

  p_death_disease_mAb114 = compute_surv(
    fit_death_disease_mAb114,
    time = state_time,
    km_limit = 14))

##pulling in induct to treat table##

par_mod <- modify(
  par_mod,

  p_tx_screen = compute_surv(
    fit_screen_to_treatment,
    time = state_time,
    km_limit = 14))

###pulling in discharge tables###

par_mod <- modify(
  par_mod,

  p_discharge_base = compute_surv(
    fit_discharge_base,
    time = state_time,
    km_limit = 14))

par_mod <- modify(
  par_mod,

  p_discharge_mAb114 = compute_surv(
    fit_discharge_mAb114,
    time = state_time,
    km_limit = 14))

par_mod <- modify(
  par_mod,

  p_discharge_REGN_EB3 = compute_surv(
    fit_discharge_REGN_EB3,
    time = state_time,
    km_limit = 14))

##This code is what imports WHO life tables##

fit_death_disease_base <- flexsurv::flexsurvreg(
  survival::Surv(time, status) ~ 1,
  dist = "weibull",
  data = tab_surv)

fit_death_disease_mAb114 <- flexsurv::flexsurvreg(
  survival::Surv(time, status) ~ 1,
  dist = "weibull",
  data = tab_surv)

fit_death_disease_REGNEB3 <- flexsurv::flexsurvreg(
  survival::Surv(time, status) ~ 1,
  dist = "weibull",
  data = tab_surv)

fit_screen_to_treatment <- flexsurv::flexsurvreg(
  survival::Surv(time, status) ~ 1,
  dist = "weibull",
  data = tab_surv)

fit_discharge_base <- flexsurv::flexsurvreg(
  survival::Surv(time, status) ~ 1,
  dist = "weibull",
  data = tab_surv)

fit_discharge_mAb114 <- flexsurv::flexsurvreg(
  survival::Surv(time, status) ~ 1,
  dist = "weibull",
  data = tab_surv)

fit_discharge_REGN_EB3 <- flexsurv::flexsurvreg(
  survival::Surv(time, status) ~ 1,
  dist = "weibull",
  data = tab_surv)

###define matrices####

##putting in point estimates from Levine et. al. paper. All diagnostic rates can be derived from these 3 data points##

sensitivity_algo = .93
specificity_algo = .23
prevalence = .42

positive_algo = sensitivity_algo * prevalence + (1- specificity_algo) * (1- prevalence)

p_pos_low = ((1-sensitivity_algo)*prevalence)/(((1- sensitivity_algo)* prevalence) + (specificity_algo * (1- prevalence)))
p_pos_hi = (sensitivity_algo * prevalence)/((sensitivity_algo * prevalence)+((1- prevalence)*(1- specificity_algo)))

##define base transition matrix##

mat_base <- define_transition(
  state_names = c("induction_discharge", "triage_hi", "triage_lo", "conf_sick", "death"),

  0, positive_algo, C, 0, 0,
  c, 0, 0, p_pos_hi*p_tx_screen, 0, 
  c, 0, 0, p_pos_low*p_tx_screen, 0, 
  p_discharge_base, 0, 0, c, p_death_disease_base,
  0, 0, 0, 0, 1)

mat_mAb_114 <- define_transition(
  state_names = c("induction_discharge", "triage_hi", "triage_lo", "conf_sick", "death"),

  0, positive_algo, C, 0, 0,
  0, 0, C, p_pos_hi*p_tx_screen, 0,
  0,0, C, p_pos_low*p_tx_screen, 0,
  p_discharge_mAb114, 0, 0, C, p_death_disease_mAb114,
  0, 0, 0, 0, 1)

mat_REGN_EB3 <- define_transition(
  state_names = c("induction_discharge", "triage_hi", "triage_lo", "conf_sick", "death"),

  0, positive_algo, C, 0, 0,
  0, 0, C, p_pos_hi*p_tx_screen, 0,
  0,0, C, p_pos_low*p_tx_screen, 0,
  p_discharge_REGN_EB3, 0, 0, C, p_death_disease_REGN_EB3,
  0, 0, 0, 0, 1)

###define state values###

##unknown represents induction - is a transient state that occurs between presentation and triage. 
##Patients do not spend a full cycle (day) there,so no qaly assigned
##some marginal ppe usage assumed as well as clinician time that always has opportunity cost in this setting

state_induction_discharge <- define_state(
  cost_total = 0, 
  qaly = 0
)

state_triage_hi <- define_state(
  cost_total = 150,
  qaly = 1/365
)

state_triage_low <- define_state(
  cost_total = 100,
  qaly = 1/365
)

Soc_cost = 200 

state_conf_sick <- define_state(
  cost_total = Soc_cost, 
  qaly = 1/365
)

###before defining the sick state, we need to allow it to only charge for mAb on day 1###

par_mod <- modify(
  par_mod,

  cost_tx_start = 300,
  cost_tx_end = Soc_cost,
  n_days = 1,
  cost_tx_cycle_mAb114 = ifelse(
    state_time < n_days,
    cost_tx_start,
    cost_tx_end)
)

state_sick_mAb114 <- define_state(
  cost_total = cost_tx_cycle_mAb114,
  qaly = 1/365
)

par_mod <- modify(
  par_mod,

  cost_tx_start = 310,
  cost_tx_end = 200,
  n_days = 1.1,
  cost_tx_cycle_REGN_EB3 = ifelse(
    state_time < n_days,
    cost_tx_start,
    cost_tx_end)
)

state_sick_REGN_Eb3 <- define_state(
  cost_total = cost_tx_cycle_REGN_EB3,
  qaly = 1/365
)

state_death <- define_state(
  cost_total = 50,
  qaly = -30
)

strat_base <- define_strategy(
  transition = mat_base,

  induction_discharge = state_induction_discharge,
  triage_hi = state_triage_hi,
  triage_lo = state_triage_low,
  conf_sick = state_conf_sick,
  death = state_death
)

strat_mAb114 <- define_strategy(
  transition = mat_mAb_114,

  induction_discharge = state_induction_discharge,
  triage_hi = state_triage_hi,
  triage_lo = state_triage_low,
  conf_sick = state_sick_mAb114,
  death = state_death
)

strat_REGN_EB3 <- define_strategy(
  transition = mat_REGN_EB3,

  induction_discharge = state_induction_discharge,
  triage_hi = state_triage_hi,
  triage_lo = state_triage_low,
  conf_sick = state_sick_REGN_Eb3,
  death = state_death
)

res_mod <- run_model(
  parameters = par_mod,

  mAb114 = strat_mAb114,
  REGN_EB3 = strat_REGN_EB3,
  base = strat_base,

  cycles = 30,

  cost = cost_total,
  effect = qaly,

  method = "life-table"
)
#> mAb114: detected use of 'state_time', expanding states: conf_sick, triage_hi, triage_lo.
#> REGN_EB3: detected use of 'state_time', expanding states: conf_sick, triage_hi, triage_lo.
#> base: detected use of 'state_time', expanding states: conf_sick, triage_hi, triage_lo.
#> Error in rep(res, nr): attempt to replicate an object of type 'builtin'

Created on 2024-01-18 with reprex v2.1.0

mmaughan0 commented 5 months ago

Understood. I pulled in excel tables for my real example but I can mock up dummy tables in the code. will get another reprex and shoot it back.

I added a generic table as a dummy in the code itself. It has the exact same error message on the original. A few issues it might be that I can think of: --the flexsurv package downloads lifetables from the WHO website. I removed the portion of the code from the paper that does this, as I am not interested in that timeframe. Its possible I removed necessary code when doing so --The matrices have time to event curves for multiple variables - its possible that at times it falls outside posssible parameters. when adding up. (i.e. cant be greater than 100% in any given cycle). I tried to set it up to not do this but not sure if I succeeded

THings I ruled out -It has nothing to do with multiplying values in the matrices. I got the same error with and without -It has nothing to do with the fact that i embedded a time to event curve in the matrices. I get the same error with and without.

Thanks so much for the help.

mmaughan0 commented 5 months ago

Update: I am getting this error (and other errors) on script I KNOW I got to run before using the heemod package. I literally went back and checked some old script line by line that I sent as progress reports to my professor. Something is going on in the background I suspect.

DavisVaughan commented 5 months ago

@mmaughan0 it looks like there are a few places in mat_base where you have a lower case c rather than an upper case C. So when it eventually tries to evaluate c, it finds the base function c(), which is a "builtin", causing the error.

If you just ensure that you are only using upper case C, it works for me

library(heemod)
library(flexsurv)
#> Loading required package: survival
library(rgho)

par_mod <- define_parameters(

    p_death_disease_base = compute_surv(
        fit_death_disease_base,
        time = state_time,
        km_limit = 14))

#establishing dummy table only for purpose of reprex. In reality each strategy would have a different table

tab_surv <- structure(list(time = c(0.4, 8.7, 7, 5.1, 9.2, 1, 0.5, 3.3, 1.8,
                                    3, 6.7, 3.7, 1.1, 5.9, 5.1, 10, 10, 10, 10, 10, 10, 10, 10, 10,
                                    10), status = c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
                                                    1L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L)), .Names = c("time",
                                                                                                                     "status"), row.names = c(NA, -25L), class = "data.frame")

##code below creates a variable in which mortality curves reside. This variable (p_death_disease_strat) 
##can then be used in transition probability matrices etc.##

par_mod <- modify(
    par_mod,

    p_death_disease_REGN_EB3 = compute_surv(
        fit_death_disease_REGNEB3,
        time = state_time,
        km_limit = 28))

par_mod <- modify(
    par_mod,

    p_death_disease_mAb114 = compute_surv(
        fit_death_disease_mAb114,
        time = state_time,
        km_limit = 14))

##pulling in induct to treat table##

par_mod <- modify(
    par_mod,

    p_tx_screen = compute_surv(
        fit_screen_to_treatment,
        time = state_time,
        km_limit = 14))

###pulling in discharge tables###

par_mod <- modify(
    par_mod,

    p_discharge_base = compute_surv(
        fit_discharge_base,
        time = state_time,
        km_limit = 14))

par_mod <- modify(
    par_mod,

    p_discharge_mAb114 = compute_surv(
        fit_discharge_mAb114,
        time = state_time,
        km_limit = 14))

par_mod <- modify(
    par_mod,

    p_discharge_REGN_EB3 = compute_surv(
        fit_discharge_REGN_EB3,
        time = state_time,
        km_limit = 14))

##This code is what imports WHO life tables##

fit_death_disease_base <- flexsurv::flexsurvreg(
    survival::Surv(time, status) ~ 1,
    dist = "weibull",
    data = tab_surv)

fit_death_disease_mAb114 <- flexsurv::flexsurvreg(
    survival::Surv(time, status) ~ 1,
    dist = "weibull",
    data = tab_surv)

fit_death_disease_REGNEB3 <- flexsurv::flexsurvreg(
    survival::Surv(time, status) ~ 1,
    dist = "weibull",
    data = tab_surv)

fit_screen_to_treatment <- flexsurv::flexsurvreg(
    survival::Surv(time, status) ~ 1,
    dist = "weibull",
    data = tab_surv)

fit_discharge_base <- flexsurv::flexsurvreg(
    survival::Surv(time, status) ~ 1,
    dist = "weibull",
    data = tab_surv)

fit_discharge_mAb114 <- flexsurv::flexsurvreg(
    survival::Surv(time, status) ~ 1,
    dist = "weibull",
    data = tab_surv)

fit_discharge_REGN_EB3 <- flexsurv::flexsurvreg(
    survival::Surv(time, status) ~ 1,
    dist = "weibull",
    data = tab_surv)

###define matrices####

##putting in point estimates from Levine et. al. paper. All diagnostic rates can be derived from these 3 data points##

sensitivity_algo = .93
specificity_algo = .23
prevalence = .42

positive_algo = sensitivity_algo * prevalence + (1- specificity_algo) * (1- prevalence)

p_pos_low = ((1-sensitivity_algo)*prevalence)/(((1- sensitivity_algo)* prevalence) + (specificity_algo * (1- prevalence)))
p_pos_hi = (sensitivity_algo * prevalence)/((sensitivity_algo * prevalence)+((1- prevalence)*(1- specificity_algo)))

##define base transition matrix##

mat_base <- define_transition(
    state_names = c("induction_discharge", "triage_hi", "triage_lo", "conf_sick", "death"),

    0, positive_algo, C, 0, 0,
    C, 0, 0, p_pos_hi*p_tx_screen, 0, 
    C, 0, 0, p_pos_low*p_tx_screen, 0, 
    p_discharge_base, 0, 0, C, p_death_disease_base,
    0, 0, 0, 0, 1)

mat_mAb_114 <- define_transition(
    state_names = c("induction_discharge", "triage_hi", "triage_lo", "conf_sick", "death"),

    0, positive_algo, C, 0, 0,
    0, 0, C, p_pos_hi*p_tx_screen, 0,
    0,0, C, p_pos_low*p_tx_screen, 0,
    p_discharge_mAb114, 0, 0, C, p_death_disease_mAb114,
    0, 0, 0, 0, 1)

mat_REGN_EB3 <- define_transition(
    state_names = c("induction_discharge", "triage_hi", "triage_lo", "conf_sick", "death"),

    0, positive_algo, C, 0, 0,
    0, 0, C, p_pos_hi*p_tx_screen, 0,
    0,0, C, p_pos_low*p_tx_screen, 0,
    p_discharge_REGN_EB3, 0, 0, C, p_death_disease_REGN_EB3,
    0, 0, 0, 0, 1)

###define state values###

##unknown represents induction - is a transient state that occurs between presentation and triage. 
##Patients do not spend a full cycle (day) there,so no qaly assigned
##some marginal ppe usage assumed as well as clinician time that always has opportunity cost in this setting

state_induction_discharge <- define_state(
    cost_total = 0, 
    qaly = 0
)

state_triage_hi <- define_state(
    cost_total = 150,
    qaly = 1/365
)

state_triage_low <- define_state(
    cost_total = 100,
    qaly = 1/365
)

Soc_cost = 200 

state_conf_sick <- define_state(
    cost_total = Soc_cost, 
    qaly = 1/365
)

###before defining the sick state, we need to allow it to only charge for mAb on day 1###

par_mod <- modify(
    par_mod,

    cost_tx_start = 300,
    cost_tx_end = Soc_cost,
    n_days = 1,
    cost_tx_cycle_mAb114 = ifelse(
        state_time < n_days,
        cost_tx_start,
        cost_tx_end)
)

state_sick_mAb114 <- define_state(
    cost_total = cost_tx_cycle_mAb114,
    qaly = 1/365
)

par_mod <- modify(
    par_mod,

    cost_tx_start = 310,
    cost_tx_end = 200,
    n_days = 1.1,
    cost_tx_cycle_REGN_EB3 = ifelse(
        state_time < n_days,
        cost_tx_start,
        cost_tx_end)
)

state_sick_REGN_Eb3 <- define_state(
    cost_total = cost_tx_cycle_REGN_EB3,
    qaly = 1/365
)

state_death <- define_state(
    cost_total = 50,
    qaly = -30
)

strat_base <- define_strategy(
    transition = mat_base,

    induction_discharge = state_induction_discharge,
    triage_hi = state_triage_hi,
    triage_lo = state_triage_low,
    conf_sick = state_conf_sick,
    death = state_death
)

strat_mAb114 <- define_strategy(
    transition = mat_mAb_114,

    induction_discharge = state_induction_discharge,
    triage_hi = state_triage_hi,
    triage_lo = state_triage_low,
    conf_sick = state_sick_mAb114,
    death = state_death
)

strat_REGN_EB3 <- define_strategy(
    transition = mat_REGN_EB3,

    induction_discharge = state_induction_discharge,
    triage_hi = state_triage_hi,
    triage_lo = state_triage_low,
    conf_sick = state_sick_REGN_Eb3,
    death = state_death
)

res_mod <- run_model(
    parameters = par_mod,

    mAb114 = strat_mAb114,
    REGN_EB3 = strat_REGN_EB3,
    base = strat_base,

    cycles = 30,

    cost = cost_total,
    effect = qaly,

    method = "life-table"
)
#> mAb114: detected use of 'state_time', expanding states: conf_sick, triage_hi, triage_lo.
#> REGN_EB3: detected use of 'state_time', expanding states: conf_sick, triage_hi, triage_lo.
#> base: detected use of 'state_time', expanding states: conf_sick, triage_hi, triage_lo.

Created on 2024-01-19 with reprex v2.0.2

mmaughan0 commented 5 months ago

Wow! I can’t believe it was something that small. Can I ask what is your process for figuring that out? Do you go follow the thread of the error message? Or try to start with a basic model and build up? I tried the latter and couldn’t get that to work either. Must have inadvertently produced the originaa add l error. Very easy to miss!!!!Sent from my iPhoneOn Jan 19, 2024, at 10:30 AM, Davis Vaughan @.***> wrote: @mmaughan0 it looks like there are a few places in mat_base where you have a lower case c rather than an upper case C. So when it eventually tries to evaluate c, it finds the base function c(), which is a "builtin", causing the error. If you just ensure that you are only using upper case C, it works for me library(heemod) library(flexsurv)

> Loading required package: survival

library(rgho)

par_mod <- define_parameters(

p_death_disease_base = compute_surv(
    fit_death_disease_base,
    time = state_time,
    km_limit = 14))

establishing dummy table only for purpose of reprex. In reality each strategy would have a different table

tab_surv <- structure(list(time = c(0.4, 8.7, 7, 5.1, 9.2, 1, 0.5, 3.3, 1.8, 3, 6.7, 3.7, 1.1, 5.9, 5.1, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10), status = c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L)), .Names = c("time", "status"), row.names = c(NA, -25L), class = "data.frame")

code below creates a variable in which mortality curves reside. This variable (p_death_disease_strat)

can then be used in transition probability matrices etc.

par_mod <- modify( par_mod,

p_death_disease_REGN_EB3 = compute_surv(
    fit_death_disease_REGNEB3,
    time = state_time,
    km_limit = 28))

par_mod <- modify( par_mod,

p_death_disease_mAb114 = compute_surv(
    fit_death_disease_mAb114,
    time = state_time,
    km_limit = 14))

pulling in induct to treat table

par_mod <- modify( par_mod,

p_tx_screen = compute_surv(
    fit_screen_to_treatment,
    time = state_time,
    km_limit = 14))

pulling in discharge tables

par_mod <- modify( par_mod,

p_discharge_base = compute_surv(
    fit_discharge_base,
    time = state_time,
    km_limit = 14))

par_mod <- modify( par_mod,

p_discharge_mAb114 = compute_surv(
    fit_discharge_mAb114,
    time = state_time,
    km_limit = 14))

par_mod <- modify( par_mod,

p_discharge_REGN_EB3 = compute_surv(
    fit_discharge_REGN_EB3,
    time = state_time,
    km_limit = 14))

This code is what imports WHO life tables

fit_death_disease_base <- flexsurv::flexsurvreg( survival::Surv(time, status) ~ 1, dist = "weibull", data = tab_surv)

fit_death_disease_mAb114 <- flexsurv::flexsurvreg( survival::Surv(time, status) ~ 1, dist = "weibull", data = tab_surv)

fit_death_disease_REGNEB3 <- flexsurv::flexsurvreg( survival::Surv(time, status) ~ 1, dist = "weibull", data = tab_surv)

fit_screen_to_treatment <- flexsurv::flexsurvreg( survival::Surv(time, status) ~ 1, dist = "weibull", data = tab_surv)

fit_discharge_base <- flexsurv::flexsurvreg( survival::Surv(time, status) ~ 1, dist = "weibull", data = tab_surv)

fit_discharge_mAb114 <- flexsurv::flexsurvreg( survival::Surv(time, status) ~ 1, dist = "weibull", data = tab_surv)

fit_discharge_REGN_EB3 <- flexsurv::flexsurvreg( survival::Surv(time, status) ~ 1, dist = "weibull", data = tab_surv)

define matrices

putting in point estimates from Levine et. al. paper. All diagnostic rates can be derived from these 3 data points

sensitivity_algo = .93 specificity_algo = .23 prevalence = .42

positive_algo = sensitivity_algo prevalence + (1- specificity_algo) (1- prevalence)

p_pos_low = ((1-sensitivity_algo)prevalence)/(((1- sensitivity_algo) prevalence) + (specificity_algo (1- prevalence))) p_pos_hi = (sensitivity_algo prevalence)/((sensitivity_algo prevalence)+((1- prevalence)(1- specificity_algo)))

define base transition matrix

mat_base <- define_transition( state_names = c("induction_discharge", "triage_hi", "triage_lo", "conf_sick", "death"),

0, positive_algo, C, 0, 0,
C, 0, 0, p_pos_hi*p_tx_screen, 0, 
C, 0, 0, p_pos_low*p_tx_screen, 0, 
p_discharge_base, 0, 0, C, p_death_disease_base,
0, 0, 0, 0, 1)

mat_mAb_114 <- define_transition( state_names = c("induction_discharge", "triage_hi", "triage_lo", "conf_sick", "death"),

0, positive_algo, C, 0, 0,
0, 0, C, p_pos_hi*p_tx_screen, 0,
0,0, C, p_pos_low*p_tx_screen, 0,
p_discharge_mAb114, 0, 0, C, p_death_disease_mAb114,
0, 0, 0, 0, 1)

mat_REGN_EB3 <- define_transition( state_names = c("induction_discharge", "triage_hi", "triage_lo", "conf_sick", "death"),

0, positive_algo, C, 0, 0,
0, 0, C, p_pos_hi*p_tx_screen, 0,
0,0, C, p_pos_low*p_tx_screen, 0,
p_discharge_REGN_EB3, 0, 0, C, p_death_disease_REGN_EB3,
0, 0, 0, 0, 1)

define state values

unknown represents induction - is a transient state that occurs between presentation and triage.

Patients do not spend a full cycle (day) there,so no qaly assigned

some marginal ppe usage assumed as well as clinician time that always has opportunity cost in this setting

state_induction_discharge <- define_state( cost_total = 0, qaly = 0 )

state_triage_hi <- define_state( cost_total = 150, qaly = 1/365 )

state_triage_low <- define_state( cost_total = 100, qaly = 1/365 )

Soc_cost = 200

state_conf_sick <- define_state( cost_total = Soc_cost, qaly = 1/365 )

before defining the sick state, we need to allow it to only charge for mAb on day 1

par_mod <- modify( par_mod,

cost_tx_start = 300,
cost_tx_end = Soc_cost,
n_days = 1,
cost_tx_cycle_mAb114 = ifelse(
    state_time < n_days,
    cost_tx_start,
    cost_tx_end)

)

state_sick_mAb114 <- define_state( cost_total = cost_tx_cycle_mAb114, qaly = 1/365 )

par_mod <- modify( par_mod,

cost_tx_start = 310,
cost_tx_end = 200,
n_days = 1.1,
cost_tx_cycle_REGN_EB3 = ifelse(
    state_time < n_days,
    cost_tx_start,
    cost_tx_end)

)

state_sick_REGN_Eb3 <- define_state( cost_total = cost_tx_cycle_REGN_EB3, qaly = 1/365 )

state_death <- define_state( cost_total = 50, qaly = -30 )

strat_base <- define_strategy( transition = mat_base,

induction_discharge = state_induction_discharge,
triage_hi = state_triage_hi,
triage_lo = state_triage_low,
conf_sick = state_conf_sick,
death = state_death

)

strat_mAb114 <- define_strategy( transition = mat_mAb_114,

induction_discharge = state_induction_discharge,
triage_hi = state_triage_hi,
triage_lo = state_triage_low,
conf_sick = state_sick_mAb114,
death = state_death

)

strat_REGN_EB3 <- define_strategy( transition = mat_REGN_EB3,

induction_discharge = state_induction_discharge,
triage_hi = state_triage_hi,
triage_lo = state_triage_low,
conf_sick = state_sick_REGN_Eb3,
death = state_death

)

res_mod <- run_model( parameters = par_mod,

mAb114 = strat_mAb114,
REGN_EB3 = strat_REGN_EB3,
base = strat_base,

cycles = 30,

cost = cost_total,
effect = qaly,

method = "life-table"

)

> mAb114: detected use of 'state_time', expanding states: conf_sick, triage_hi, triage_lo.

> REGN_EB3: detected use of 'state_time', expanding states: conf_sick, triage_hi, triage_lo.

> base: detected use of 'state_time', expanding states: conf_sick, triage_hi, triage_lo.

Created on 2024-01-19 with reprex v2.0.2

—Reply to this email directly, view it on GitHub, or unsubscribe.You are receiving this because you were mentioned.Message ID: @.***>

DavisVaughan commented 5 months ago

I ran your reprex in RStudio, used traceback() to see where the error came from:

> traceback()
7: FUN(X[[i]], ...)
6: lapply(x_tidy, function(x) {
       res <- rlang::eval_tidy(x, data = p2)
       if (length(res) == 1) {
           return(rep(res, nr))
       }
       res
   })
5: eval_transition.uneval_matrix(uneval_transition, parameters)
4: eval_transition(uneval_transition, parameters)
3: eval_strategy(strategy = uneval_strategy_list[[n]], parameters = parameters, 
       init = init, cycles = cycles, method = method, expand_limit = state_time_limit[[n]], 
       inflow = inflow, strategy_name = n)
2: run_model_(uneval_strategy_list = uneval_strategy_list, parameters = parameters, 
       init = init, cycles = cycles, method = method, cost = new_quosure(substitute(cost), 
           env = parent.frame()), effect = new_quosure(substitute(effect), 
           env = parent.frame()), state_time_limit = state_time_limit, 
       central_strategy = central_strategy, inflow = inflow)
1: run_model(parameters = par_mod, mAb114 = strat_mAb114, REGN_EB3 = strat_REGN_EB3, 
       base = strat_base, cycles = 30, cost = cost_total, effect = qaly, 
       method = "life-table")

Then did debugonce(heemod:::eval_transition.uneval_matrix) and ran it again and then stepped through that function until it got to the following line

6: lapply(x_tidy, function(x) {
       res <- rlang::eval_tidy(x, data = p2)
       if (length(res) == 1) {
           return(rep(res, nr))
       }
       res
   })

Then I investigated that until I figured out which iteration was actually failing, and I looked at it in x_tidy and saw it had the lower case c in it.

(Note that I actually had to call debugonce(heemod:::eval_transition.uneval_matrix) a few times, because that function gets called 3 times here and it is only on the 3rd one that we see the failure)

mmaughan0 commented 5 months ago

That traceback function I did not know. Amazing. Question - do you take tips? Assuming you do this for free, but you figured out way quicker than anyone I could hire (I tried). Saved me probably a week of futile attempts at debugging. Sent from my iPhoneOn Jan 19, 2024, at 10:52 AM, Davis Vaughan @.***> wrote: I ran your reprex in RStudio, used traceback() to see where the error came from:

traceback() 7: FUN(X[[i]], ...) 6: lapply(x_tidy, function(x) { res <- rlang::eval_tidy(x, data = p2) if (length(res) == 1) { return(rep(res, nr)) } res }) 5: eval_transition.uneval_matrix(uneval_transition, parameters) 4: eval_transition(uneval_transition, parameters) 3: eval_strategy(strategy = uneval_strategy_list[[n]], parameters = parameters, init = init, cycles = cycles, method = method, expand_limit = state_time_limit[[n]], inflow = inflow, strategy_name = n) 2: runmodel(uneval_strategy_list = uneval_strategy_list, parameters = parameters, init = init, cycles = cycles, method = method, cost = new_quosure(substitute(cost), env = parent.frame()), effect = new_quosure(substitute(effect), env = parent.frame()), state_time_limit = state_time_limit, central_strategy = central_strategy, inflow = inflow) 1: run_model(parameters = par_mod, mAb114 = strat_mAb114, REGN_EB3 = strat_REGN_EB3, base = strat_base, cycles = 30, cost = cost_total, effect = qaly, method = "life-table") Then did debugonce(heemod:::eval_transition.uneval_matrix) and ran it again and then stepped through that function until it got to the following line 6: lapply(x_tidy, function(x) { res <- rlang::eval_tidy(x, data = p2) if (length(res) == 1) { return(rep(res, nr)) } res }) Then I investigated that until I figured out which iteration was actually failing, and I looked at it in x_tidy and saw it had the lower case c in it. (Note that I actually had to call debugonce(heemod:::eval_transition.uneval_matrix) a few times, because that function gets called 3 times here and it is only on the 3rd one that we see the failure)

—Reply to this email directly, view it on GitHub, or unsubscribe.You are receiving this because you were mentioned.Message ID: @.***>

DavisVaughan commented 5 months ago

Assuming you do this for free

We work for Posit (https://posit.co/), so we do get paid to do this! The best thing you can do is to keep using our free software, and to tell your company / school / friends about our paid products (https://posit.co/products/enterprise/)! 😄