You signed in with another tab or window. Reload to refresh your session.You signed out in another tab or window. Reload to refresh your session.You switched accounts on another tab or window. Reload to refresh your session.Dismiss alert
You can consult reprex below - issue is that DSA will not run with multiple variables. I have tried many variables individually but the second I add another one I get this error:
#> Error in [.data.frame(x$complete_parameters, 1, dsa$variables): undefined columns selected
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: 0x000001c03f2b3c80>#> <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 statereprex#> 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: 0x000001c04b442c98>#> <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 tablestime=state_time, ##state time means time in that given state, in this case it will be in the "conf-sick" statekm_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 treatmentpar_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 probabilitiesfit_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=.93specificity_algo=.23prevalence=.42## the below just rearranges all the terms and definitions above algebraically. no new information is here, just rearrangingpositive_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 estimate0, 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 cost0, 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 estimate0, 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 cost0, 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 estimate0, 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 cost0, 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 inductionstate_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 locationsPer_Survivor_Cost=830Per_Dead_Cost=321Avg_2014_Pt_Cost= (Per_Survivor_Cost/11.2+Per_Dead_Cost/4.2)/2CPI_Multiplier_2014_2024=1.32Avg_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 triagePer_Survivor_Cost_low=446Per_Dead_Cost_low=185Avg_2014_Pt_Cost_low= (Per_Survivor_Cost_low/11.2+Per_Dead_Cost_low/4.2)/2Avg_2024_SoC_Pt_Cost_low=CPI_Multiplier_2014_2024*Avg_2014_Pt_Cost_lowCI_Per_Surv_cost=862-800CI_Per_Dead_cost=351-292CI_Per_Surv_cost_low=466-428CI_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=3702BF_CN_ann_wage=2823Zambia_RN_ann_wage=6476Nigeria_RN_ann_wage=8709Ghana_RN_ann_wage=7680BF_ratio_RN_CN=BF_RN_ann_wage/BF_CN_ann_wageZambia_CN_ann_wage=Zambia_RN_ann_wage/BF_ratio_RN_CNNigeria_CN_ann_wage=Nigeria_RN_ann_wage/BF_ratio_RN_CNGhana_CN_ann_wage=Ghana_RN_ann_wage/BF_ratio_RN_CN##Step 4: do weighted average of CN/RN wage based on wilson paperRN_pct_DC=78/102Ghana_avg_wage=RN_pct_DC*Ghana_RN_ann_wage+ (1-RN_pct_DC)*Ghana_CN_ann_wageZambia_avg_wage=RN_pct_DC*Zambia_RN_ann_wage+ (1-RN_pct_DC)*Zambia_CN_ann_wageNigeria_avg_wage=RN_pct_DC*Nigeria_RN_ann_wage+ (1-RN_pct_DC)*Nigeria_CN_ann_wageBF_avg_wage=RN_pct_DC*BF_RN_ann_wage+ (1-RN_pct_DC)*BF_CN_ann_wageRegional_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 paperTotal_ETU_Patients=600Confirmed_cases=prevalence*Total_ETU_PatientsTraining_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/191Mild_DW=.004Moderate_Prev=2/191Moderate_DW=.033Blindness_Prev=10/191Blindness_DW=.195DW_Vision_Impair=Mild_Prev*Mild_DW+Moderate_Prev*Moderate_DW+Blindness_Prev*Blindness_DWstate_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 inputTLC_Cost=6890Days_Operation=120Resupply_Cost_ETU=92Amortization_Rate=Days_Operation/(365*5)
Refrigeration_Cost=TLC_Cost*Amortization_RateCCL_cost_ETU=Resupply_Cost_ETU+Refrigeration_CostCCL_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 pointcost_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 mAb114def_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
You can consult reprex below - issue is that DSA will not run with multiple variables. I have tried many variables individually but the second I add another one I get this error:
#> Error in
[.data.frame(x$complete_parameters, 1, dsa$variables): undefined columns selectedNote 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)
Created on 2024-03-06 with reprex v2.1.0
The text was updated successfully, but these errors were encountered: