── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr 1.1.2 ✔ readr 2.1.4
✔ forcats 1.0.0 ✔ stringr 1.5.0
✔ ggplot2 3.4.3 ✔ tibble 3.2.1
✔ lubridate 1.9.2 ✔ tidyr 1.3.0
✔ purrr 1.0.2
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag() masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(heemod)
Attaching package: 'heemod'
The following object is masked from 'package:purrr':
modify
library(diagram)
Loading required package: shape
library(ggplot2)
Build parameter for base model
## import the Canada 2020 life table with mortality value in all ages (both sex)life_table<-read.csv('./life_table_CAN_2020.csv')# remove row 111 from life_table (age older than 111 )life_table <- life_table %>%filter(Age.group !="110 years and over") %>%separate(Age.group, into =c("Age", "Unit"), sep =" ") %>%select(Age, VALUE) %>%mutate(Age=as.numeric(Age), fatality_daily =rescale_prob(p = VALUE, from =365 ))# build age parameter,par_allerg <-define_parameters(age_initial =1,age =floor(age_initial + model_time/365))# age-dependent remission probabilityar_age<-function(age) {if_else(age <=6, rescale_prob(0.058,from =365), 0) }#_ww: watch and wait scenario#_ED: transfer to ED scenario #ns: non severe#sw: severe watch and wait,#sED: severe transfer to ED,#sh: severe with hospitalized, #ar: food allergy remission, #faf: food allergy fatality,#dr: discount rate#acm: all-cause # parameter for X10 mortality in watch and wait #parameter for transition #parameter for transition par_allerg<-modify(par_allerg,p_severe =0.087, # transition from non-severe reaction to severe reaction p_ns_ar_original =0.058,p_ns_ar =if_else(age <=6, rescale_prob(p_ns_ar_original,from =365), 0), # transition from non-severe reaction to food allergy remissionp_biphas =0.046, # have biphasitic reactionp_EAI =0.906, # patient with known food allergy having an EAIp_ns_sw_ww = (p_EAI-p_biphas) *rescale_prob(p =p_severe, from =365) , # non-severe reaction to severe reaction for watch and waitp_ns_sED_ww = (1-p_EAI+p_biphas) *rescale_prob(p = p_severe, from =365),# non-severe reaction to severe reaction transfer to ED for watch and waitp_faf =0.00000069, # food allergy fatality among patients with food allergyp_sh_faf =0.0045, # transition from hospitalization to food allergy fatalityp_sh_ED =rescale_prob(p_faf, from =365)/rescale_prob(p_severe, from =365)/p_sh_faf, # transition from severe state ED to hospitalization p_ns_sED_ED =rescale_prob(p =p_severe, from =365), # transition from non-severe to ED in ED scenarioacm =look_up(data = life_table, Age = age, #daily all-cause mortalityvalue ="fatality_daily") )# parameter for cost and utilitypar_allerg<-modify(par_allerg,p_EAI_ED =0.54, # Probability of using an EAI prior to ED visit treatment_cost_ED =0.8+95*p_EAI_ED, # epi treatment cost in ED treatment_cost_ww =95, # EAI treatment cost in watch and wait treatment_cost_ww_ED =0.8, # epi cost for ED transfer in watch and wait scenario remission_cost_all =569/365, # remission patient daily cost for all scenariomedical_cost_ns_original =1393, # yearly direct medical cost of food allergy < 18 yrs medical_cost_original =304, # Average cost of all-cause ED visitmedical_cost_ED_ED =round(medical_cost_original *1.09,digits =0), # daily medical cost in ED no_visit =0.3, # no.of ED visit per yearmedical_cost_ns =round(((medical_cost_ns_original*1.068)-(medical_cost_ED_ED *no_visit )), digits =0)/365,# daily medical cost for Canadian with food allergyambulance_cost_ED_ED =848, # ambulance cost medical_h =3670, # total medical cost hospitalization for allergylength_stay =2.1, # length of stay per admission medical_cost_hospital =round(1.068*medical_h/length_stay,digits =0), # medical cost hospitalizedindirect_cost_allergy_original =4173, # indirect cost for food allergy < 18 yearsout_of_pocket_allergy_original =2440, # out of pocket cost of food allergy per year indirect_cost_hospital =251, # indirect cost in hospitalization - societal out_of_pocket_original =89, # out of pocket cost out_of_pocket_emergency =round(out_of_pocket_original *1.068, digits =0), # out of pocket in ED-societal waiting_time_ED =3.6, # waiting time in ED/ hours indirect_cost_emergency =round(waiting_time_ED*31.37,digits =0), # indirect cost in ED -societal indirect_cost_allergy =round(1.068*(indirect_cost_allergy_original-no_visit* indirect_cost_emergency), digits =0) /365, # indirect cost for food allergy-societalout_of_pocket_cost_allergy =round(1.068* (out_of_pocket_allergy_original-no_visit* indirect_cost_emergency), digits =0) /365, # out of pocket cost for food allergy-societal utilify_far_original =0.93,utility_far = utilify_far_original /365, # utility -food allergy remission utility_fa_original =0.92,utility_sa =0.91,utility_fa = utility_fa_original/365, # utility- food allergy non-severeutility_sr = utility_fa -(1-utility_sa)/365# utility-severe reaction )
## Model X 100 mortality#parameter for #parameter for X100 mortality in watch and wait par_allerg_100<-modify( par_allerg,p_ww_faf = (100*rescale_prob(p_faf, from =365)/rescale_prob(p_severe, from =365) -(1-p_EAI+p_biphas)* p_sh_ED*p_sh_faf)/(p_EAI-p_biphas),# Watch and wait to food allergy fatality X100, baseline 0.000002dr =rescale_prob( p=0.015, from =365) )allergy_mod_100<-run_model(parameters = par_allerg_100,ED_transfer = strategy_ED,watch_wait = strategy_watch,init = time0,cycles =20*365,cost = cost_total,effect = utility_total)summary(allergy_mod_100)
## Model X 500 mortality#parameter for #parameter for X500 mortality in watch and wait par_allerg_500<-modify( par_allerg,p_ww_faf = (500*rescale_prob(p_faf, from =365)/rescale_prob(p_severe, from =365) -(1-p_EAI+p_biphas)* p_sh_ED*p_sh_faf)/(p_EAI-p_biphas), # Watch and wait to food allergy fatality X500, baseline 0.000002dr =rescale_prob( p=0.015, from =365) )allergy_mod_500<-run_model(parameters = par_allerg_500,ED_transfer = strategy_ED,watch_wait = strategy_watch,init = time0,cycles =20*365,cost = cost_total,effect = utility_total)summary(allergy_mod_500)
## Model X 1000 mortality# parameter for X 1000 mortality in watch and waitpar_allerg_1000<-modify( par_allerg,p_ww_faf = (1000*rescale_prob(p_faf, from =365)/rescale_prob(p_severe, from =365) -(1-p_EAI+p_biphas)* p_sh_ED*p_sh_faf)/(p_EAI-p_biphas), # Watch and wait to food allergy fatality X500, baseline 0.000002dr =rescale_prob( p=0.015, from =365) )allergy_mod_1000<-run_model(parameters = par_allerg_1000,ED_transfer = strategy_ED,watch_wait = strategy_watch,init = time0,cycles =20*365,cost = cost_total,effect = utility_total)summary(allergy_mod_1000)
# run base model for discount = 0% par_allerg_dr_0<-modify( par_allerg,p_ww_faf = (10*rescale_prob(p_faf, from =365)/rescale_prob(p_severe, from =365) -(1-p_EAI+p_biphas)* p_sh_ED*p_sh_faf)/(p_EAI-p_biphas), #transition from watch and wait to food allergy fatality dr =0)allergy_mod_dr_0<-run_model(parameters = par_allerg_dr_0,ED_transfer = strategy_ED,watch_wait = strategy_watch,init = time0,cycles =20*365,cost = cost_total,effect = utility_total)summary(allergy_mod_dr_0)
Get output result and table after simulation (with out PSA)
value_10<-get_values(allergy_mod_10)counts_10<-get_counts(allergy_mod_10)#develop the table to compare the effect of different mortality in watch & wait#function to build tablevaluetable<-function(value,counts, senario_name){ a <-value %>%group_by(.strategy_names, value_names) %>%summarize(.groups ="keep", sum=sum(value/10000)) %>%pivot_wider(names_from = value_names, values_from = sum) %>%ungroup() %>%mutate(cost_diff = cost_total-lag(cost_total)) %>%mutate(Qaly_diff = utility_total-lag(utility_total)) %>%mutate(cost_per_QALY_saved= cost_diff/Qaly_diff) %>%mutate(scenario =senario_name)b <-counts %>%filter(state_names =="state_faf"& model_time ==7300) %>%mutate(death_diff = (count -lag(count))/10000)merge(a,b, by =".strategy_names" ) %>%mutate(cost_per_life_saved =abs(cost_diff/death_diff)) %>%select(scenario, .strategy_names, cost_total,utility_total,cost_diff,Qaly_diff,cost_per_QALY_saved, cost_per_life_saved,death_diff ) }value_table_10<-valuetable(value_10,counts_10, "X10_mortality_in_ww")#INMB calculationINMB<-value_table_10 %>%mutate(INMB_result = Qaly_diff*50000-cost_diff )print(INMB)
scenario .strategy_names cost_total utility_total cost_diff
1 X10_mortality_in_ww ED_transfer 21797.71 15.90335 NA
2 X10_mortality_in_ww watch_wait 20640.86 15.90262 -1156.85
Qaly_diff cost_per_QALY_saved cost_per_life_saved death_diff
1 NA NA NA NA
2 -0.0007276453 1589854 12586613 9.191115e-05
INMB_result
1 NA
2 1120.468
#Build the final table to compare the cost_per_lifesaved with different mortalityfinal_Table<-value_table_10 %>%mutate_if(is.numeric, as.character) final_Table <-replace(final_Table, is.na(final_Table),"reference") final_Table<-final_Table %>%rename("strategy_name"=".strategy_names","QALYs_total"="utility_total","increment_cost"="cost_diff","increment_Qaly"="Qaly_diff","increment_cost_per_life_saved"="cost_per_life_saved")#final summary table for different fatality scenarioprint(as.data.frame(final_Table))
