mmaughan0 commented 5 months ago

You can consult reprex below - issue is that DSA will not run with multiple variables. I have tried many variables individually and the code runs just fine but the second I add another one I get this error:

> Error in [.data.frame(x$complete_parameters, 1, dsa$variables): undefined columns selected

I consulted Flipovic-Perucci, Zarca and Isabelle Durand-Zaleski's Paper "Markov Models for Health Ecomonc Evaluations: THe R Package heemod" for correct syntax.

Note that many of these values are dummies that I threw in just to get the reprex to run without errors (excluding the final one of course)

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: 0x000002175d4ecc30>
#> <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: 0x000002176bbccce8>
#> <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,

  n_days = 11,
  death_transition_prob_base = ifelse(
    state_time < n_days, 
    p_death_disease_base, 
    0)
)

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, 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(
  Avg_2024_SoC_Pt_Cost, 0, 10000,
  Tx_cost, 0, 15000)

res_dsa <- run_dsa(res_mod, dsa = def_dsa)
#> Running DSA on strategy 'mAb114'...
#> Error in `[.data.frame`(x$complete_parameters, 1, dsa$variables): undefined columns selected

^{Created on 2024-03-06 with reprex v2.1.0}

KZARCA commented 5 months ago

Ok, it's a frequent error: all your variables must be defined within the function define_parameters, and not in the global environment

mmaughan0 commented 5 months ago

ok - I ran into that issue before with the PSA so I understand what you mean - but in that DSA both of those variables work individually. Is it really the same issue?

KZARCA commented 5 months ago

yes I modified your code locally, and it's working

mmaughan0 commented 5 months ago

ok - do you have any insight as to why it would run one at a time?

mmaughan0 commented 5 months ago

Ok now I am confusing myself a little bit. So we have both PSA and DSA, and you are saying that both need to be within a modify parameter statement, however

PSA Corrections and comment: 1) I corrected my issue with the PSA by taking a variable OUT of the par_mod <- modify statement, specifically 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) ) -

The PSA would not work if Tx_Cost was modified in the statement above. It did work when I added Tx_Cost to the statement below, which as you will notice is NOT in a parmod statement state_sick_REGN_Eb3 <- define_state( cost_total = Tx_cost + CCL_Cost_pp + Training_cost_pp + Avg_2024_SoC_Pt_Cost, qaly = 0)

DSA Corrections and Comment: I have Tx_Cost and SoC Cost within a par_mod statement and not within a par_mod statement, however it only runs when I have it in the statement below one at a time. The error I am getting does not seem to have anything to do with how I define the parameter.

def_dsa <- define_dsa( Tx_cost, 0, 10000)

Thanks!

KZARCA commented 5 months ago

The issue is with these variables: 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_DCGhana_RN_ann_wage + (1-RN_pct_DC)Ghana_CN_ann_wage Zambia_avg_wage = RN_pct_DCZambia_RN_ann_wage + (1-RN_pct_DC)Zambia_CN_ann_wage Nigeria_avg_wage = RN_pct_DCNigeria_RN_ann_wage + (1-RN_pct_DC)Nigeria_CN_ann_wage BF_avg_wage = RN_pct_DCBF_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

Regional_daily_wage = Regional_avg_wage/working_days_per_year

Total_ETU_Patients = 600 Confirmed_cases = prevalence* Total_ETU_Patients

Training_cost_pp = Regional_daily_wage*102/Confirmed_cases

LY_point_est = 43.2

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_PrevMild_DW + Moderate_PrevModerate_DW + Blindness_Prev*Blindness_DW

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/(3655) Refrigeration_Cost = TLC_Cost Amortization_Rate CCL_cost_ETU = Resupply_Cost_ETU + Refrigeration_Cost

CCL_Cost_pp = CCL_cost_ETU/Confirmed_cases

KZARCA commented 5 months ago

which are defined outside define_parameters

mmaughan0 commented 5 months ago

ok - I am going to dig into this and put all of these variables into a par_mod statement - will report back if I cannot get it to work.

Would it be fair to say that the fact that only one variable would work at a time in the DSA is just a downstream effect of not definining the variables within a parameter? If I fix the problem it shouldnt matter, but I am just curious at this point. It seems an interesting side effect.

KZARCA commented 5 months ago

It's not a side effect, really. That is the way heemod is coded: relying heavily on non standard evaluation and capturing the calling environment.

aphp / heemod

DSA will not run multiple variables in R heemod #12

> Error in [.data.frame(x$complete_parameters, 1, dsa$variables): undefined columns selected

we are going to do same calculation for the not as severe supportive care from Bartsch 2014 so we can estimate SoC low triage

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

Step 4: do weighted average of CN/RN wage based on wilson paper

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

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