KZARCA
commented
6 months ago Did you use death_transition_prob_base_sense somewhere in your model?
Closed mmaughan0 closed 6 months ago
KZARCA
commented
6 months ago Did you use death_transition_prob_base_sense somewhere in your model?
mmaughan0
commented
6 months ago I did on the original but forgot to add into the reprex
library(heemod)
library(flexsurv)
#> Loading required package: survival
library(rgho)
library(dplyr)
#>
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#>
#> filter, lag
#> The following objects are masked from 'package:base':
#>
#> intersect, setdiff, setequal, union
library(readxl)
library(reprex)
##setting up directory, as we have many tables to put in that are too large to practically code##
setwd("C:/Users/Matt/OneDrive/Documents/Portfolio Project")
##reading in necessary tables from directory##
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: 0x0000020c948e0330>
#> <environment: namespace:base>
##par_mod <- define_parameters(
##age_base = 30*365,##cycle in this model is 1 day
##age_cycle = model_time + age_base)
## I muted this code above because I don't think its necessary to import WHO life tables - as my cycle length is 1 day and the lifetables
#are based on years I do not want to muddle the two. I can realize the lost life years during transition to the death state
reprex
#> function (x = NULL, input = NULL, wd = NULL, venue = c("gh",
#> "r", "rtf", "html", "slack", "so", "ds"), render = TRUE,
#> advertise = NULL, session_info = opt(FALSE), style = opt(FALSE),
#> comment = opt("#>"), tidyverse_quiet = opt(TRUE), std_out_err = opt(FALSE),
#> html_preview = opt(TRUE), outfile = deprecated(), show = deprecated(),
#> si = deprecated())
#> {
#> if (lifecycle::is_present(show)) {
#> html_preview <- show
#> lifecycle::deprecate_warn(when = "1.0.0", what = "reprex(show)",
#> with = "reprex(html_preview)")
#> }
#> if (lifecycle::is_present(si)) {
#> session_info <- si
#> }
#> reprex_impl(x_expr = substitute(x), input = input, wd = wd,
#> venue = venue, render = render, new_session = TRUE, advertise = advertise,
#> session_info = session_info, style = style, html_preview = html_preview,
#> comment = comment, tidyverse_quiet = tidyverse_quiet,
#> std_out_err = std_out_err, outfile = outfile)
#> }
#> <bytecode: 0x0000020c995cb6e8>
#> <environment: namespace:reprex>
par_mod <- define_parameters(
p_death_disease_base = compute_surv(
fit_death_disease_base, ##fit death disease base will be defined later, and will pull in TTE tables
time = state_time, ##state time means time in that given state, in this case it will be in the "conf-sick" state
km_limit = 28)) ## km_limit is 28 because the tables I have go to 28 days
##The above code is calculating probability based on a life curve that I will define in a later piece of code##
par_mod <- modify(
par_mod,
p_death_disease_mAb114 = compute_surv(
fit_death_disease_mAb114,
time = state_time,
km_limit = 28))
##the above code does evertything the p_death_disease_base does, but for the REGNEB3 treatment
par_mod <- modify(
par_mod,
p_death_disease_REGN_EB3 = compute_surv(
fit_death_disease_REGNEB3,
time = state_time,
km_limit = 28))
## the above code does everything p_death_disease_base does, but for the mAb114 treatment
##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))
## the above code is setting a probability of moving to the next stage for the time between induction(screen) and treatment.
##fit_screen_treatment will be defined later
###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_REGNEB3 = compute_surv(
fit_discharge_REGNEB3,
time = state_time,
km_limit = 28))
## the above probabilities function just like the death probabilities, but for discharge
##dummy table to allow the reprex to work. actual model has different tables for each##
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")
##The below code is what fits the curves into probabilities
fit_death_disease_base <- flexsurv::flexsurvreg(
survival::Surv(time, status) ~ 1,
dist = "gompertz",
data = tab_surv)
fit_death_disease_mAb114 <- flexsurv::flexsurvreg(
survival::Surv(time, status) ~ 1,
dist = "gompertz",
data = tab_surv)
fit_death_disease_REGNEB3 <- flexsurv::flexsurvreg(
survival::Surv(time, status) ~ 1,
dist = "gompertz",
data = tab_surv)
fit_screen_to_treatment <- flexsurv::flexsurvreg(
survival::Surv(time, status) ~ 1,
dist = "gompertz",
data = tab_surv)
fit_discharge_base <- flexsurv::flexsurvreg(
survival::Surv(time, status) ~ 1,
dist = "gompertz",
data = tab_surv)
fit_discharge_mAb114 <- flexsurv::flexsurvreg(
survival::Surv(time, status) ~ 1,
dist = "gompertz",
data = tab_surv)
fit_discharge_REGNEB3 <- flexsurv::flexsurvreg(
survival::Surv(time, status) ~ 1,
dist = "gompertz",
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##
# prevalence may correlate with sens/spec in an epidemiological sense, but not in a mathematical sense. For our purposes they are considered independent
#Definitions:
#TP = True Positive
#TN = True Negative
#FP = False Positive
#FN = False Negative
#Sensitivity = TP/(TP+FN) --> Think "If I am positive, what are the chances I test positive"
#Specificity = TN/(TN+FP) --> Think "If I am negative, what are the chances I test negative"
# Prevalence = (TP+FN)/(TP+TN+FP+FN) --> What are the chances a person in the given population has the disease
##drew the below from Levine et. al. 2015 - which helpfully has confidence intervals, and thus can be coded probabilistically
sensitivity_algo = .93
specificity_algo = .23
prevalence = .42
## the below just rearranges all the terms and definitions above algebraically. no new information is here, just rearranging
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##
par_mod <- modify(
par_mod,
screen_delay = 0
)
par_mod <- modify(
par_mod,
n_days = 11,
death_transition_prob_base = ifelse(
state_time < n_days,
p_death_disease_base,
0)
)
par_mod <- modify(
par_mod,
death_transition_prob_base_sense = ifelse(
screen_delay == 0,
death_transition_prob_base,
death_transition_prob_base * (1.11)^screen_delay)
)
par_mod <- modify(
par_mod,
n_days = 11,
death_transition_prob_mAb114 = ifelse(
state_time < n_days,
p_death_disease_mAb114,
0)
)
par_mod <- modify(
par_mod,
n_days = 11,
death_transition_prob_REGNEB3 = ifelse(
state_time < n_days,
p_death_disease_REGN_EB3,
0)
)
par_mod <- modify(
par_mod,
n_days = 11,
screen_transition_prob = ifelse(
state_time < n_days - 5,
p_tx_screen,
1)
)
mat_base <- define_transition(
state_names = c("induction", "triage_hi", "triage_lo", "conf_sick", "death", "terminal"),
0, positive_algo, C, 0, 0, 0, ##you can only go from induction (where you start) to two states - hi or low triage. the positive algo gives a point estimate
0, C, 0, screen_transition_prob*p_pos_hi, 0, screen_transition_prob*(1-p_pos_hi), ##p_tx_screen*p_pos_hionce in triage_hi, the p_tx_screen gives the chances of moving to next stage by day, but multiplies this by chance of being positive.
0, 0, C, screen_transition_prob*p_pos_low, 0, screen_transition_prob*(1-p_pos_low), ##p_tx_screen*p_pos_low same as above but starting from low
0, 0, 0, C, death_transition_prob_base_sense, p_discharge_base, ##there are two curves here from the confirmed sick state, you can either die or be discharged. the "plug" is remaining in that state. I am honestly not sure if this will work.
0, 0, 0, 0, 0, 1, ##death covers death through burial. it is not a terminal state - leaving it as a terminal state would count the death costs every remaining cycle, vastly overcounting QALY loss and burial cost
0, 0, 0, 0, 0, 1) #terminal is the only terminal state.
plot(fit_death_disease_mAb114)
mat_mAb_114 <- define_transition(
state_names = c("induction", "triage_hi", "triage_lo", "conf_sick", "death", "terminal"),
0, positive_algo, C, 0, 0, 0, ##you can only go from induction (where you start) to two states - hi or low triage. the positive algo gives a point estimate
0, C, 0, screen_transition_prob*p_pos_hi, 0, screen_transition_prob*(1-p_pos_hi), ##p_tx_screen*p_pos_hionce in triage_hi, the p_tx_screen gives the chances of moving to next stage by day, but multiplies this by chance of being positive.
0, 0, C, screen_transition_prob*p_pos_low, 0, screen_transition_prob*(1-p_pos_low), ##p_tx_screen*p_pos_low same as above but starting from low
0, 0, 0, C, death_transition_prob_mAb114, p_discharge_mAb114, ##there are two curves here from the confirmed sick state, you can either die or be discharged. the "plug" is remaining in that state. I am honestly not sure if this will work.
0, 0, 0, 0, 0, 1, ##death covers death through burial. it is not a terminal state - leaving it as a terminal state would count the death costs every remaining cycle, vastly overcounting QALY loss and burial cost
0, 0, 0, 0, 0, 1) #terminal is the only terminal state.
mat_REGN_EB3 <- define_transition(
state_names = c("induction", "triage_hi", "triage_lo", "conf_sick", "death", "terminal"),
0, positive_algo, C, 0, 0, 0, ##you can only go from induction (where you start) to two states - hi or low triage. the positive algo gives a point estimate
0, C, 0, screen_transition_prob*p_pos_hi, 0, screen_transition_prob*(1-p_pos_hi), ##p_tx_screen*p_pos_hionce in triage_hi, the p_tx_screen gives the chances of moving to next stage by day, but multiplies this by chance of being positive.
0, 0, C, screen_transition_prob*p_pos_low, 0, screen_transition_prob*(1-p_pos_low), ##p_tx_screen*p_pos_low same as above but starting from low
0, 0, 0, C, death_transition_prob_REGNEB3, p_discharge_REGNEB3, ##there are two curves here from the confirmed sick state, you can either die or be discharged. the "plug" is remaining in that state. I am honestly not sure if this will work.
0, 0, 0, 0, 0, 1, ##death covers death through burial. it is not a terminal state - leaving it as a terminal state would count the death costs every remaining cycle, vastly overcounting QALY loss and burial cost
0, 0, 0, 0, 0, 1) #terminal is the only terminal state.
###define state values###
state_induction <- define_state(
cost_total = 0,
qaly = 0
)
## costs are still dummy values, BUT the following should be included
##SoC costs (Bartsh et. al) (all strategies, all states within the ETU)
##incremental cold chain costs (Mvundura et. al. WHO/UNICEF) (experimental arms)
## Incremental Training costs for mAb (conf_sick state in experimental arms, first day, spread across all pos patients)
##burial costs (for death state, all arms)
## qalys are negligible in every state but death, where the lost life years over the course of a lifetime are considered
## can absolutely code probabilistically based on life tables, have found literature of average age at induction
state_triage_hi <- define_state(
cost_total = Avg_2024_SoC_Pt_Cost,
qaly = 0
)
state_triage_low <- define_state(
cost_total = Avg_2024_SoC_Pt_Cost_low,
qaly = 0
)
state_terminal <- define_state(
cost_total = 0,
qaly = 0
)
## step 1 - Bartsch et. al. 2014 defines ppe and drug cost per episode for both those who recover and those who die, and duration of treatment.
##step 2: divide this cost for those who recovered (830/11.2) = 74.1 and those who die 321/4.2 = 76.4
##step 3: average out to 75
##step 4: convert to 2024 dollars.https://www.bls.gov/data/inflation_calculator.htm Decided to use exchange rates rather than PPP due to Bartsch methodology and variety of possible payors and locations
Per_Survivor_Cost = 830
Per_Dead_Cost = 321
Avg_2014_Pt_Cost = (Per_Survivor_Cost/11.2 + Per_Dead_Cost/4.2)/2
CPI_Multiplier_2014_2024 = 1.32
Avg_2024_SoC_Pt_Cost = CPI_Multiplier_2014_2024 * Avg_2014_Pt_Cost
### we are going to do same calculation for the not as severe supportive care from Bartsch 2014 so we can estimate SoC low triage
Per_Survivor_Cost_low = 446
Per_Dead_Cost_low = 185
Avg_2014_Pt_Cost_low = (Per_Survivor_Cost_low/11.2 + Per_Dead_Cost_low/4.2)/2
Avg_2024_SoC_Pt_Cost_low = CPI_Multiplier_2014_2024 * Avg_2014_Pt_Cost_low
CI_Per_Surv_cost = 862-800
CI_Per_Dead_cost = 351-292
CI_Per_Surv_cost_low = 466-428
CI_Per_Dead_cost_low = 202-169
##step 1: Wilson 2014 suggests ~100 patients on a continuing basis (80 when she arrived, 120 at peak). She mentions 102 direct care personnel (78/24 registered vs. certified) total
##step 2: assuming incremental training on mAbs will take all of these nurses 1 day of wages
##step 3: average wage calc - McCoy et al. 2008 - 4 different countries - assuming registered nurse closer to the 78 nurses and certified closer to assistants
##3a: Burkina Faso 2454 RN 1871 CN (2006) Ghana RN 4896 (2005) Zambia 4128 (2005) Nigeria 5097 (2001). Use BF ratio for other countries
## the below have already been run through the CPI calculator in Jan of their respective year to Jan 23##
BF_RN_ann_wage = 3702
BF_CN_ann_wage = 2823
Zambia_RN_ann_wage = 6476
Nigeria_RN_ann_wage = 8709
Ghana_RN_ann_wage = 7680
BF_ratio_RN_CN = BF_RN_ann_wage/BF_CN_ann_wage
Zambia_CN_ann_wage = Zambia_RN_ann_wage/BF_ratio_RN_CN
Nigeria_CN_ann_wage = Nigeria_RN_ann_wage/BF_ratio_RN_CN
Ghana_CN_ann_wage = Ghana_RN_ann_wage/BF_ratio_RN_CN
##Step 4: do weighted average of CN/RN wage based on wilson paper
RN_pct_DC = 78/102
Ghana_avg_wage = RN_pct_DC*Ghana_RN_ann_wage + (1-RN_pct_DC)*Ghana_CN_ann_wage
Zambia_avg_wage = RN_pct_DC*Zambia_RN_ann_wage + (1-RN_pct_DC)*Zambia_CN_ann_wage
Nigeria_avg_wage = RN_pct_DC*Nigeria_RN_ann_wage + (1-RN_pct_DC)*Nigeria_CN_ann_wage
BF_avg_wage = RN_pct_DC*BF_RN_ann_wage + (1-RN_pct_DC)*BF_CN_ann_wage
Regional_avg_wage = mean(Ghana_avg_wage, Zambia_avg_wage, Nigeria_avg_wage, BF_avg_wage)
working_days_per_year = 250
##5 days per week, 50 weeks per year, accounting for travel time sickness etc.##
Regional_daily_wage = Regional_avg_wage/working_days_per_year
##step 5: add wage back up for all health workers, allocate it at a patient level, using total patients in WIlson paper
Total_ETU_Patients = 600
Confirmed_cases = prevalence* Total_ETU_Patients
Training_cost_pp = Regional_daily_wage*102/Confirmed_cases
###DEFINE STATES####
#First define DALYs so that we can place into distribution for sensitivity analysis later#
LY_point_est = 43.2
class(LY_dist)
#> Error in eval(expr, envir, enclos): object 'LY_dist' not found
###All Prev Figures taken from Shantha et. al.###
###All DW Figures taken from Salomon et. al.###
Mild_Prev = 14/191
Mild_DW = .004
Moderate_Prev = 2/191
Moderate_DW = .033
Blindness_Prev = 10/191
Blindness_DW = .195
DW_Vision_Impair = Mild_Prev*Mild_DW + Moderate_Prev*Moderate_DW + Blindness_Prev*Blindness_DW
state_sick_base <- define_state(
cost_total = Avg_2024_SoC_Pt_Cost,
qaly = 0
)
###before defining the sick state, we need to allow it to only charge for mAb on day 1###
##going to first define Soc per day cost##
Tx_cost = 500
##this is a discreet decision point less than a probabilistic sensitivity. the types of people using this model as a decision tool often are also negotiating the price, so it can be an output of the decsion making process, not an input
TLC_Cost= 6890
Days_Operation = 120
Resupply_Cost_ETU = 92
Amortization_Rate = Days_Operation/(365*5)
Refrigeration_Cost = TLC_Cost * Amortization_Rate
CCL_cost_ETU = Resupply_Cost_ETU + Refrigeration_Cost
CCL_Cost_pp = CCL_cost_ETU/Confirmed_cases
##McCoy et. al.
par_mod <- modify(
par_mod,
cost_tx_start = Avg_2024_SoC_Pt_Cost + Tx_cost + CCL_Cost_pp + Training_cost_pp,
cost_tx_end = Avg_2024_SoC_Pt_Cost,
n_days = 1,
cost_tx_cycle_mAb114 = ifelse(
state_time < n_days-9.9,
cost_tx_start,
cost_tx_end)
)
state_sick_mAb114 <- define_state(
cost_total = cost_tx_cycle_mAb114,
qaly = 0
)
par_mod <- modify(
par_mod,
cost_tx_start = Avg_2024_SoC_Pt_Cost + Tx_cost + CCL_Cost_pp + Training_cost_pp,
cost_tx_end = Avg_2024_SoC_Pt_Cost,
n_days = 11,
cost_tx_cycle_REGN_EB3 = ifelse(
state_time < n_days-9.9, ## we first defined in n days in the probability matrix to cut off the KM curve at 10 days. For that reason we have to use it as a reference point
cost_tx_start,
cost_tx_end)
)
state_sick_REGN_Eb3 <- define_state(
cost_total = cost_tx_cycle_REGN_EB3,
qaly = 0
)
state_death <- define_state(
cost_total = 50,
qaly = (LY_point_est*-1)*(1-DW_Vision_Impair)
)
##define strategies##
strat_base <- define_strategy(
transition = mat_base,
induction = state_induction,
triage_hi = state_triage_hi,
triage_lo = state_triage_low,
conf_sick = state_sick_base,
death = state_death,
terminal = state_terminal
)
strat_mAb114 <- define_strategy(
transition = mat_mAb_114,
induction = state_induction,
triage_hi = state_triage_hi,
triage_lo = state_triage_low,
conf_sick = state_sick_mAb114,
death = state_death,
terminal = state_terminal
)
strat_REGN_EB3 <- define_strategy(
transition = mat_REGN_EB3,
induction = state_induction,
triage_hi = state_triage_hi,
triage_lo = state_triage_low,
conf_sick = state_sick_REGN_Eb3,
death = state_death,
terminal = state_terminal
)
res_mod <- run_model(
parameters = par_mod,
mAb114 = strat_mAb114,
REGN_EB3 = strat_REGN_EB3,
base = strat_base,
cycles = 40,
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.
summary(res_mod, threshold = c(1000,5000,15000))
#> 3 strategies run for 40 cycles.
#>
#> Initial state counts:
#>
#> induction = 1000L
#> triage_hi = 0L
#> triage_lo = 0L
#> conf_sick = 0L
#> death = 0L
#> terminal = 0L
#>
#> Counting method: 'life-table'.
#>
#> Values:
#>
#> cost_total qaly
#> mAb114 935086.6 -7664.1
#> REGN_EB3 935086.6 -7664.1
#> base 719247.6 -7664.1
#>
#> Net monetary benefit difference:
#>
#> 1000 5000 15000
#> mAb114 0.000 0.000 0.000
#> REGN_EB3 0.000 0.000 0.000
#> base 215.839 215.839 215.839
#>
#> Efficiency frontier:
#>
#> base
#>
#> Differences:
#>
#> Cost Diff. Effect Diff. ICER Ref.
#> REGN_EB3 0.000 0 NaN mAb114
#> base -215.839 0 -Inf mAb114
def_dsa <- define_dsa(
screen_delay, 0, 5)
res_dsa <- run_dsa(res_mod, dsa = def_dsa)
#> Running DSA on strategy 'mAb114'...
#> Running DSA on strategy 'REGN_EB3'...
#> Running DSA on strategy 'base'...
plot(res_dsa, type = "simple", result = "effect", strategy= "base")
Created on 2024-03-07 with reprex v2.1.0
mmaughan0
commented
6 months ago now it worked in the reprex but it doesnt work in my original.... back to teh drawing board I guess
mmaughan0
commented
6 months ago Ok this is very strange - I figured out the issue but I am not sure how to fix it. If you leave the base case for screen_delay at zero - it will not vary in the DSA - however, if base case is above zero it WILL vary. the problem is I want my base case to be zero.
mmaughan0
commented
6 months ago Found a work around - base case is very small number close to zero, and the conditional statement is directional
par_mod <- modify( par_mod,
screen_delay = .1
)
par_mod <- modify( par_mod,
n_days = 11, death_transition_prob_base = ifelse( state_time < n_days, p_death_disease_base, 0) )
par_mod <- modify( par_mod,
death_transition_prob_base_sense = ifelse( screen_delay <=.1, death_transition_prob_base, death_transition_prob_base * (1.11)^screen_delay)
I want to embed an indicator within the transition probability matrix that allows mortality to vary with a change in screening delay. I can get the point estimate to change and the DSA to run - but the DSA does not return any variation in effect when I do so. I wrapped everything in the define_parameters function which was previous error I had. FUll reprex is below but I have the most relevant code immediately below
par_mod <- modify( par_mod,
screen_delay = 0
)
par_mod <- modify( par_mod,
death_transition_prob_base_sense = ifelse( screen_delay == 0, death_transition_prob_base, death_transition_prob_base * (1.11)^screen_delay) )
Created on 2024-03-07 with reprex v2.1.0