Markov cost effectiveness model - survival curve is using a best fit line - I do not want it to. #1909
Comments
|
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 |
|
Cmd/Ctrl + V? |
|
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 foundCreated on 2024-01-17 with reprex v2.1.0 |
|
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 |
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!!! |
|
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! |
|
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. |
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 |
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: THings I ruled out Thanks so much for the help. |
|
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. |
|
@mmaughan0 it looks like there are a few places in If you just ensure that you are only using upper case 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 |
|
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: ***@***.***>
|
|
I ran your reprex in RStudio, used > 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 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 (Note that I actually had to call |
|
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: 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)
—Reply to this email directly, view it on GitHub, or unsubscribe.You are receiving this because you were mentioned.Message ID: ***@***.***>
|
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/)! 😄 |
mmaughan0
changed the title
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
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"
)
The text was updated successfully, but these errors were encountered: