--- title: "Cost-Effectiveness of delaying procession in AD: Aducanumab, Lecanemab, and Donaenemab" author: "Minseon Chung" date: today output: pdf_document --- ```{r set up} knitr::opts_chunk$set(echo = TRUE) library(tidyverse) library(ggplot2) library(knitr) library(gt) library(rmarkdown) library(quantreg) library(boot) library(stats) library(rigr) library(dplyr) library(tidyr) library(grDevices) library(dampack) library(kableExtra) ``` ###Basecase Analysis### ##Defining Parameters #1. Model Settings# ```{r} n_cohort <- 1 #simulation cohort size time_horizon_yr <- 10 #time_horizon in years time_horizon_mnths <- 10*12 weekspermonth <- 4.333 weeksperyear <- 52 #number of weeks per year time_horizon_wks <- ceiling(time_horizon_mnths*weekspermonth) cycle <- 1 #cycle length = 1 week cycle_unit <- "weeks" n_cycles <- (time_horizon_yr*weeksperyear)+1 #weeks #Set sex variable sex_2cat <- c("Male","Female") n_sex <- length(sex_2cat) #For DSA param_range <- 0.2 #20% #Transition States s_names <- c("MCI","Mild AD","Moderate AD","Severe AD","Dead") n_states <- length(s_names) s_names_4cat <- c("MCI","Mild AD","Moderate AD","Severe AD") n_states_4cat <- length(s_names_4cat) #Treatments t_names <- c("SOC","Aducanumab", "Lecanemab", "Donanemab") #SOC=standard of care n_strategies <- length(t_names) #Setting setting_names <- c("Community","Institutionalized") n_setting <- length(setting_names) #Discounts dc_costs <- 0.03 dc_outcomes <- 0.03 ``` #2. Inputs# ##A. Baseline Population ```{r} base_age <- 71.1 # base_male_perc <- 0.46 #male percentage base_female_perc <- 1-base_male_perc base_cdrsb <- 3.3 ##A. Health States (PMID: 18695059) cdrsb_raw <- base_cdrsb cdrsb <- floor(cdrsb_raw*2)/2 #cdrsb_raw is rounded to nearest 0.5. # ad_stages <- cut(cdrsb, # breaks = c(0,0.5,4.5,9.5,16.0,18.0), # labels = c("Normal","MCI","Mild AD","Moderate AD","Severe AD"), # right = FALSE, #use [left,right) intervals # include.lowest = TRUE) ad_stage <- function(cdrsb_val){ stage <- case_when( cdrsb_val >= 0 & cdrsb_val < 0.5 ~ "Normal", cdrsb_val >= 0.5 & cdrsb_val < 4.5 ~ "MCI", cdrsb_val >= 4.5 & cdrsb_val < 9.5 ~ "Mild AD", cdrsb_val >= 9.5 & cdrsb_val < 16.0 ~ "Moderate AD", cdrsb_val >= 16.0 ~ "Severe AD") return(stage) } ``` ##B. Costs ```{r} #Percentage of population > 65yo cover under medicaid medicaid_perc <- if(cdrsb < 4.0) 0.08 else if(cdrsb >= 4.0 & cdrsb <= 5.0) 0.16 else if(cdrsb > 5.5 & cdrsb <= 7.0) 0.25 else if(cdrsb > 7.5 & cdrsb <= 9.0) 0.34 else if(cdrsb > 9.5 & cdrsb <= 12.0) 0.39 else if(cdrsb > 12.0) 0.62 stage <- c(NA) #Patient Healthcare per month (in USD $) cost_pt_healthcare <- c( "MCI" = 1291.83, "Mild AD" = 1515.16, "Moderate AD" = 2010.99, "Severe AD" = 2309.11 ) #Caregiver(care partner) Healthcare per month (in USD $) cost_cp_healthcare <- c( "MCI" = 788.32, "Mild AD" = 817.48, "Moderate AD" = 834.14, "Severe AD" = 846.64 ) #Patient Non-healthcare per month (in USD $) cost_pt_nonhealthcare <- c( "MCI" = 233.52, "Mild AD" = 431.44, "Moderate AD" = 684.85, "Severe AD" = 1149.80 ) #Patient Out of Pocket Cost per month (in USD $) #pmid 37863109 page 51. patient oop may be different by residential settings? cost_pt_oop <- c( "MCI" = 176.18, "Mild AD" = 176.18, "Moderate AD" = 176.18, "Severe AD" = 176.18 ) #Patient Productivity per month (in USD $) cost_pt_productivity <- c( "MCI" = 2628.19, "Mild AD" = 1624.96, "Moderate AD" = 920.20, "Severe AD" = 0 ) #Caregiver Lost Productivity per month (in USD $): Community & institutionalized (monthly hours*hourly wage) us_hourly_wage <- 34.75 #May 2025 cost_cp_lostprod_community <- c( "MCI" = 69 * us_hourly_wage, "Mild AD" = 113 * us_hourly_wage, "Moderate AD" = 126 * us_hourly_wage, "Severe AD" = 298 * us_hourly_wage ) cost_cp_lostprod_institute <- c( "MCI" = 30 * us_hourly_wage, "Mild AD" = 50 * us_hourly_wage, "Moderate AD" = 74 * us_hourly_wage, "Severe AD" = 131 * us_hourly_wage ) #Long Term Care Cost (per month) cost_ltc <- c( "MCI" = 10461.44, "Mild AD" = 10461.44, "Moderate AD" = 10461.44, "Severe AD" = 10461.44 ) #$9,698 per month according to pmid 38331952 (2020) #https://www.carescout.com/cost-of-care #9,277 according to administration on aging ``` ##C. Institutionalization ```{r} ltc_percent <- c( "MCI" = 0.017, "Mild AD" = 0.043, "Moderate AD" = 0.116, "Severe AD" = 0.432 ) ``` ##D. Utilities ```{r} # Utility (Community Setting) utility_comm <- c( "MCI" = 0.73, "Mild AD" = 0.68, "Moderate AD" = 0.54, "Severe AD" = 0.37 ) # Utility (Long-Term Care Setting) utility_ltc <- c( "MCI" = 0.73, "Mild AD" = 0.71, "Moderate AD" = 0.48, "Severe AD" = 0.31 ) # Care Partner Disutility disutility_cp <- c( "MCI" = 0.03, "Mild AD" = 0.05, "Moderate AD" = 0.08, "Severe AD" = 0.10 ) # Care Partner Utility utility_cp <- c( "MCI" = 0.88, "Mild AD" = 0.76, "Moderate AD" = 0.65, "Severe AD" = 0.49 ) ``` ##E. Mortality*** ```{r} #AD realted mortality mortality_ad_value = c( "MCI" = 1.04, #assumed 1, because not significantly different from the referenced CN. "Mild AD" = 2.36, "Moderate AD" = 2.75, "Severe AD" = 8.04 ) #PMID: 25697699, 37580727(2023)(used Exploratory HR (adjusted with all predictors and covariates); Curta used 9921853 (1999) #All-Cause Background Mortality (Human Mortality Database http://mortality.org/) background_mortality <- read.csv("C:/Users/MinC/Desktop/UW/Project AD Progression/R/background_mortality.csv") mortality_combined_perwk <- function(ad_stage, age, sex){ mort_AD_HR <- if(ad_stage == "MCI"){ mortality_ad_value["MCI"] }else if(ad_stage == "Mild AD"){ mortality_ad_value["Mild AD"] }else if(ad_stage == "Moderate AD"){ mortality_ad_value["Moderate AD"] }else if(ad_stage == "Severe AD"){ mortality_ad_value["Severe AD"] }else(stop("Error")) a <- floor(age) #background_mortality[row,column] mort_background <- if(a < 100 & sex == "Male"){ (background_mortality[a+1,4]) } else if(a < 100 & sex == "Female"){ (background_mortality[a+1,5]) } else if(a >= 100 & sex == "Male"){ (background_mortality[101,4]) } else if(a >= 100 & sex == "Female"){ (background_mortality[101,5]) } else{ stop("Enter appropriate age and sex") } #convert annual probability of death to weekly probability mort_rate <- -log(1-mort_background) mort_prob_pwk <- 1- exp(-mort_rate/weeksperyear) mort_probability <- 1-(1-mort_prob_pwk)^mort_AD_HR mort_probability <- mort_AD_HR*mort_background mort_probability <- mort_probability/weeksperyear return(mort_probability) } ``` ##F. Treatment Specific Parameters ```{r} #Clinical Efficacy# #Aducanumab endpoint_t_adu <- 18#primary endpoint time or trial length tx_efficacy_adu <- 0.22 #primary endpoint reduction in progression #tx_eff_stop_adu <- #cdr-sb level at which treatment benefit stops effect_adu <- tx_efficacy_adu/endpoint_t_adu #Lecanemab endpoint_t_leca <- 18#primary endpoint time tx_efficacy_leca <- 0.2711 #primary endpoint reduction in progression effect_leca <- tx_efficacy_leca/endpoint_t_leca #Donanemab endpoint_t_dona <- 17.5#primary endpoint time tx_efficacy_dona <- 0.289 #primary endpoint reduction in progression effect_dona <- tx_efficacy_dona/endpoint_t_dona #Drug Cost Per Month (USD)# cost_txtitration_pmonth_adu <- 479.40 #per 1.7ml bottle. multiply x1.76 for month 3,4(3ml bottle), and x3.53 for month 5,6 (3ml bottle). Then only use cost_tx_pmonth_adu afterward cost_tx_pmonth_adu <- 2171.40 #300mg/3ml #85mg for 1,2 + 255mg for 3,4 + 510mg for 5,6 iv cost_tx_pmonth_leca <- 2327.78 #=2*(646.61+(2*258.64)) RedBook #does not require titration cost_tx_pmonth_dona <- 2782.60 #4*695.65 RedBook #titration: 700mg for iv 1,2,3. #Administration Cost Per Month (USD)# CPT 96365 + 96366 cost_admin_pmonth_adu <- 57.90 #Aducanumab, iv over an ~hour [PI]: 96366*3, 2nd+3rd: 96366*2, 4th: 96366 cost_admin_pmonth_leca <- 57.90*2 #Lecanemab, iv over an ~hour [PI]: 1st - 96366*3, 2nd+3rd: 96366*2, 4th: 96366 cost_admin_pmonth_dona <- 57.90 #Donanemab, iv over 30 min [PI]: 1st-4th iv - 96366, 5th, no 96366. #Donanemab appropriate use recommendations suggest up to 1hr post-infusion monitoring for first 4 infusions and 30 min post-iv monitoring the following infusions. #Lecanemab AUR suggest 3 hrs after the 1st, 2hrs after the 2nd and 3rd, and 30 min for following infusions. ``` #G. Treatment Discontinuation ```{r} #Discontinuation Rate (per month) # dc_rate_pmonth_adu <- 0.0914/endpoint_t_adu #9.14% d/c'ed adu during the 18 month trial # dc_rate_pmonth_leca <- 0.069/endpoint_t_leca #6.9% d/c'ed leca during the 18 month trial # dc_rate_pmonth_dona <- 0.131/endpoint_t_dona #13.1% d/c'ed dona during the 17.5 month trial # #Discontinuation divided into 6 months. And applied only for the first 6 months of starting tx.(ICER 2022 Beta-amyloid CE) dc_rate_6month_adu <- 0.0914/6 #9.14% d/c'ed adu during the 18 month trial dc_rate_6month_leca <- 0.069/6 #6.9% d/c'ed leca during the 18 month trial dc_rate_6month_dona <- 0.131/6 #13.1% d/c'ed dona during the 17.5 month trial # #Discontinuation divided into 12 months. And applied only for the 1st year of starting tx.(ICER 2022 Beta-amyloid CE) # dc_rate_1YR_adu <- 0.0914/12 #9.14% d/c'ed adu during the 18 month trial # dc_rate_1YR_leca <- 0.069/12 #6.9% d/c'ed leca during the 18 month trial # dc_rate_1YR_dona <- 0.131/12 #13.1% d/c'ed dona during the 17.5 month trial ``` #H. Monitoring ```{r} mri_cost_70553 <- 314.08 #retreived 05/04/25 monitor_adu <- function(week, treatment){ if(treatment == 1){ monitor_cost <- 0 }else if(treatment == 2 & week %in% c(0, 20, 28, 36, 48)){ #PMID: 34585212, Package insert #Aducanumab: at week 0, 20, 28, 36, 48 (prior to 5th, 7th, 9th and 12th infusion; high risK: + 10th dose) monitor_cost <- mri_cost_70553 }else if(treatment == 3 & week %in% c(0, 10, 14, 28)){ #PMID: 37357276 #Lecanemab: at week 0, 10, 14, 28 (after 0th, 5th, 7th, 14th infusion; high risk: + week 52) monitor_cost <- mri_cost_70553 }else if(treatment == 4 & week %in% c(0, 8, 12, 16, 28)){ #PMID: 40155270 #Donanemab: at week 0, 8, 12, 16, 28 (prior to 2nd, 3rd, 4th, 7th infusion; high risk: 12th) monitor_cost <- mri_cost_70553 }else(monitor_cost <- 0) return(monitor_cost) } ``` #I. Adverse Events ```{r} #Assumption on ADE: Serious ADE only occurs in the first 6 month of starting tx's. There is chance the providers would d/c meds after ADE for first few months. #Management Cost +/- HCPCS in USD per event drvisit_cost_99213 <- 88.95 #CMS 2025 drug_migraine_cost <- 19.22 #PMID 30142988 mri12wks_cost <- 314.08*3 #every month for 12 weeks add_postiv_monitoring_96366 <- 20.31 #$20.31 = CPT 96366 Additional 30 min monitoring following initial post-infusion monitoring (CMS 2024B ) #Symptomatic ARIA (moderate ARIA typically resolve in 12 wks) aria_pweek_adu <- 0.2444/(6*weekspermonth) #ex. 24.44% of pt had symptomatic ARIA at 18 month aria_pweek_leca <- 0.1606/(6*weekspermonth) aria_pweek_dona <- 0.1656/(6*weekspermonth) aria_disutility <- 0.14 #distuility from headache. Persist for 7 days #won't account for aria disutilities, since it's already taken account into headache disutility. aria_mngm_cost <- drvisit_cost_99213*3 + mri12wks_cost #doctor visit + mri every month for 12 weeks #Mild headache headache_pweek_adu <- 0.2009/(6*weekspermonth) headache_pweek_leca <- 0.1114/(6*weekspermonth) headache_pweek_dona <- 0.1395/(6*weekspermonth) headache_disutility <- 0.14 #PMID: 30711061 (SuppAppex pg . HA lasts 2 days) headache_mngm_cost <- drug_migraine_cost #Fall fall_pweek_adu <- 0.1465/(6*weekspermonth) fall_pweek_leca <- 0.1036/(6*weekspermonth) fall_pweek_dona <- 0.1337/(6*weekspermonth) fall_disutility <- 0.2 #Persist for 7 days fall_mngm_cost <- drvisit_cost_99213 #Infusion Related Reactions #typically resolves in 24hours irr_pweek_adu <- 0.0298/(6*weekspermonth) irr_pweek_leca <- 0.2639/(6*weekspermonth) irr_pweek_dona <- 0.0868/(6*weekspermonth) irr_disutility <- 0.011*1 #PMID: 20224930; Mild Injection Reaction -0.011 (SD: 0.158) #PMID 30522373 #disutility = 0.066 per IRR event. IRR last for 24hours. irr_mngm_cost <- 20.31 #$20.31 = CPT 96366 Additional 30 min monitoring following initial post-infusion monitoring #Diphenhydramine 50mg or acetaminophen. Very miniscule amount a few dollars (see below) #In case of IRR: The diphenhydramine or acetaminophen may be repeated every 4–6 hours until symptoms fully resolve.(PMID: 37357276) #Diphenhydramine: 25 to 50 mg; after 15 to 30 minutes, if symptoms persist, may repeat dose as needed. Do not exceed 100 mg within a 1-hour period. Potential dose: 50 mg PO Q6HR = 200mg -> $0.1596 (Micromedex) #Acetaminophen: for infusion-related pain/fever): 325 to 650 mg every 4 to 6 hours as needed or 1 g every 6 hours as needed; maximum dose: 4 g/day; Potential dose: 500mg PO Q6HT = 2000mg -> $0.0096 (Micromedex) ####Disutility related to 1-2 hour infusion # disu_iv <- 0.04??? ``` #J. Perspective Inclusion ```{r} # commerical <- 1 #1 to include, 0 to exclude # medicare <- 1 #1 to include, 0 to exclude # medicaid <- 1 #1 to include, 0 to exclude # ltc <- 1 #1 to include, 0 to exclude#long-term care # oop <- 1 #1 to include, 0 to exclude#out of pocket # pt_productivity <- 1 #1 to include, 0 to exclude # caregiver_productivity <- 1 #1 to include, 0 to exclude # caregiver_healthcare <- 1 #1 to include, 0 to exclude ``` #K.1. Transition Probabilities: Kuhnel et al 2021 (PMID: 33853199) ```{r} ########################Kuhnel et al 2021 pmid: 33853199########################## #Standard of Care l_soc <- 0.1030 g_soc <- 3.6900 v_soc <- 0.2750 #Define time at 0? t_0 <- log((base_cdrsb-v_soc)/l_soc)*exp(g_soc) #Aducanumab delta_adu_soc <- (exp((t_0+endpoint_t_adu)/exp(g_soc))-exp(t_0/exp(g_soc)))*l_soc g_adu <- log(endpoint_t_adu/log(delta_adu_soc*((1-tx_efficacy_adu)/(base_cdrsb-v_soc))+1)) l_adu <- (base_cdrsb-v_soc)/exp(t_0/exp(g_adu)) v_adu <- 0.275 #Lecanemab delta_leca_soc <- (exp((t_0+endpoint_t_leca)/exp(g_soc))-exp(t_0/exp(g_soc)))*l_soc g_leca <- log(endpoint_t_leca/log(delta_leca_soc*((1-tx_efficacy_leca)/(base_cdrsb-v_soc))+1)) l_leca <- (base_cdrsb-v_soc)/exp(t_0/exp(g_leca)) v_leca <- 0.275 #Donanemab delta_dona_soc <- (exp((t_0+endpoint_t_dona)/exp(g_soc))-exp(t_0/exp(g_soc)))*l_soc g_dona <- log(endpoint_t_dona/log(delta_dona_soc*((1-tx_efficacy_dona)/(base_cdrsb-v_soc))+1)) l_dona <- (base_cdrsb-v_soc)/exp(t_0/exp(g_dona)) v_dona <- 0.275 kuhnel_param <- data.frame( treatment = c(1,2,3,4), #1: SOC, 2: Aducanumab, 3: Donanemab, 4: Lecanemab l = c(l_soc, l_adu, l_leca, l_dona), g = c(g_soc, g_adu, g_leca, g_dona), v = c(v_soc, v_adu, v_leca, v_dona) ) #Setting cycles months <- 0:((n_cycles-1)/weekspermonth) weeks <- 0:(n_cycles-1) years <- 0:(((n_cycles-1)/weekspermonth)/12) #function to calculate CDR-SB at time t kuhnel_cdrsb <- function(set_cycle, treatment, kuhnel_param){ # Sele0ct the correct parameters for the given treatment params <- kuhnel_param[kuhnel_param$treatment == treatment, ] l <- params$l g <- params$g v <- params$v cycle_name <- deparse(substitute(set_cycle)) t_0 <- log((base_cdrsb-v_soc)/l_soc)*exp(g_soc) # Convert t_cycle into month unit for the equation if (cycle_name == "weeks") { t <- (set_cycle/4.333)+t_0 } else if (cycle_name == "months") { t <- set_cycle + t_0 } else if (cycle_name == "years") { t <- (set_cycle*12) + (t_0) } else { stop("Invalid time unit. Choose 'weeks', 'months', or 'years'.") } # Predicted CDR-SB with tx benefit stopping once pts reach "Moderate AD" (CDR-SB = 9.5) predicted_cdrsb <- (l * exp(t / exp(g)) + v) soc_cdrsb <- (l_soc * exp(t/ exp(g_soc)) + v_soc) if(treatment == 1){ # Set max cdr-sb (18.0) and convert them to 18.0 max_cdrsb <- which(predicted_cdrsb >= 18.0) predicted_cdrsb[max_cdrsb] <- 18.0 # Retrieve result cdrsb <- predicted_cdrsb }else{ #for treatment = 2,3, and 4 # Convert the predicted_cdrsb trajectory by adding no month benefit effect after CDR-SB > 9.5 + no more progression after CDR-SB = 18.0 # 1. Find the month where no_benefit starts valid_indices <- which(predicted_cdrsb > 9.5) time_nobenefit <- valid_indices[1] new_time_nobenefit <-if (cycle_name == "weeks") { (time_nobenefit/4.333) + t_0 } else if (cycle_name == "months") { time_nobenefit + t_0 } else if (cycle_name == "years") { (time_nobenefit*12) + t_0 } else { stop("Error on month_nobenefit") } # 2. Apply the finding from (1) and convert predict_cdrsb accordingly nobenefit_cdrsb <- which(predicted_cdrsb >= 9.5) cdrsb_delta <- ((l_soc * exp(new_time_nobenefit/ exp(g_soc)) + v_soc) - (l * exp(new_time_nobenefit/ exp(g)) + v)) predicted_cdrsb[nobenefit_cdrsb] <- soc_cdrsb[nobenefit_cdrsb] - cdrsb_delta #SOC_cdrsb - CDR_delta # nobenefit_cdrsb <- which(predicted_cdrsb < 18.0 & predicted_cdrsb >= 9.5) # cdrsb_delta <- ((l_soc * exp(new_t/ exp(g_soc)) + v_soc) - (l * exp(new_t/ exp(g)) + v)) # predicted_cdrsb[nobenefit_cdrsb] <- (l_soc * exp(new_t / exp(g_soc)) + v_soc) - cdrsb_delta # #SOC_cdrsb - CDR_delta # Set max cdr-sb (18.0) and convert them to 18.0 max_cdrsb <- which(predicted_cdrsb >= 18.0) predicted_cdrsb[max_cdrsb] <- 18.0 # if(predicted_cdrsb >= 18.0){ # predicted_cdrsb <- 18.0 # }else if(predicted_cdrsb >= 9.5){ # # Retrieve the month where the tx benefit stops # valid_indices <- which(predicted_cdrsb > 9.5) # month_nobenefit <- valid_indices[1] # # # Retrieve CDRSB_delta @ month no benefit # cdr_delta <- (l_soc * exp(t_0 + month_nobenefit/ exp(g_soc)) + v_soc) - (l * exp(t_0 + month_nobenefit/ exp(g)) + v) # predicted_cdrsb <- (l_soc * exp(t / exp(g_soc)) + v_soc) - cdr_delta # }else{ # predicted_cdrsb <- predicted_cdrsb # } # Retrieve Result cdrsb <- predicted_cdrsb } return(cdrsb) } #l_soc*exp((t_0)/exp(g_soc))+v_soc #l_adu*exp((t_0)/exp(g_adu))+v_adu cdrsb_soc <- kuhnel_cdrsb(weeks, 1, kuhnel_param) cdrsb_adu <- kuhnel_cdrsb(weeks, 2, kuhnel_param) cdrsb_leca <- kuhnel_cdrsb(weeks, 3, kuhnel_param) cdrsb_dona <- kuhnel_cdrsb(weeks, 4, kuhnel_param) plot(weeks, cdrsb_soc, type = "l", lwd = 2, xlab = "Weeks", ylab = "CDR-SB", main = "CDR-SB Progression for SOC", col = "grey", ylim = c(0,18), yaxt = "n") + lines(weeks, cdrsb_adu, lwd = 2, col = "red") + lines(weeks, cdrsb_leca, lwd = 2, col = "blue") + lines(weeks, cdrsb_dona, lwd = 2, col = "green") + axis(side = 2, at = seq(0,18,by =1)) #plot(months, cdrsb_adu, type = "l", lwd = 2, xlab = "Months", ylab = "CDR-SB", # main = "CDR-SB Progression for Aducanumab") #plot(months, cdrsb_leca, type = "l", lwd = 2, xlab = "Months", ylab = "CDR-SB", # main = "CDR-SB Progression for Lecanemab") #plot(months, cdrsb_dona, type = "l", lwd = 2, xlab = "Months", ylab = "CDR-SB", # main = "CDR-SB Progression for Donanemab") write.csv(cdrsb_soc, "cdrsb_soc1.csv", row.names = FALSE) write.csv(cdrsb_adu, "cdrsb_adu1.csv", row.names = FALSE) write.csv(cdrsb_leca, "cdrsb_leca1.csv", row.names = FALSE) write.csv(cdrsb_dona, "cdrsb_dona1.csv", row.names = FALSE) ``` #K.2 Transitional Trajectory for Tx D/C'ed Patients mid-therapy ```{r} #kuhnel_cdrsb_dc predicts the progression trajectory that pts w/ tx discontinuation follows, depending on the time they d/c their treatment (time_dc) #treatment - only applied to tx 2,3 and 4. #j serves as at the time of d/c in the loop #set_cycle should be same everywhere = "Weeks" in this analysis #treatment in numerics = i, time_dc = time of dc kuhnel_cdrsb_dc <- function(treatment, time_dc){ # Redefine treatment for convenience i <- treatment # Call the equation params appropriate for the treatment params <- kuhnel_param[kuhnel_param$treatment == i, ] l <- params$l g <- params$g v <- params$v # Convert t_cycle into month unit for the d/c trajectory t_0 <- log((base_cdrsb-v_soc)/l_soc)*exp(g_soc) t <- (weeks/4.333)+t_0 # Retrieve reference trajectory old_traj <-(l * exp(t/ exp(g)) + v) # Predicted CDR-SB with tx benefit stopping once pts reach "Moderate AD" (CDR-SB = 9.5) soc_cdrsb <- (l_soc * exp(t/ exp(g_soc)) + v_soc) # Define cycles after tx d/c cycle_postdc <- time_dc:521 # Convert the predicted_cdrsb trajectory by adding no month benefit effect after CDR-SB > 9.5 + no more progression after CDR-SB = 18.0 time_nobenefit <- cycle_postdc[1] # 2. Apply the finding from (1) and convert predict_cdrsb accordingly cdrsb_delta <- soc_cdrsb[time_nobenefit] - old_traj[time_nobenefit] old_traj[cycle_postdc] <- soc_cdrsb[cycle_postdc] - cdrsb_delta #SOC_cdrsb - CDR_delta # Set max cdr-sb (18.0) and convert them to 18.0 max_cdrsb <- which(old_traj >= 18.0) old_traj[max_cdrsb] <- 18.0 # Set max cdr-sb (18.0) and convert them to 18.0 max_cdrsb <- which(soc_cdrsb >= 18.0) soc_cdrsb[max_cdrsb] <- 18.0 cdrsb <- if(i %in% c(2,3,4)){old_traj}else if(i == 1){soc_cdrsb} return(cdrsb) } cdrsb_4_dcatwk5 <- kuhnel_cdrsb_dc(4,300) cdrsb_2_dcatwk5 <-kuhnel_cdrsb_dc(2,300) plot(weeks, cdrsb_4_dcatwk5, type = "l", lwd = 2, xlab = "Weeks", ylab = "CDR-SB", main = "CDR-SB Progression for SOC", col = "grey", ylim = c(0,18), yaxt = "n") + lines(weeks, cdrsb_2_dcatwk5, lwd = 2, col = "red") ``` ############################ MARKOV ANALYSIS ############################ #1. Plot Progression Trajectories ```{r} ####Plot Progression Trajectories##### #CDR-SB Progression Curves by treatment (1 = SOC, 2 = Aducanumab, 3 = Lecanemab, 4 = Donanemab) #No clinical benefit applied after pts enters "Moderate AD" (CDR-SB >= 9.5) cdrsb_soc <- kuhnel_cdrsb(weeks, 1, kuhnel_param) cdrsb_adu <- kuhnel_cdrsb(weeks, 2, kuhnel_param) cdrsb_leca <- kuhnel_cdrsb(weeks, 3, kuhnel_param) cdrsb_dona <- kuhnel_cdrsb(weeks, 4, kuhnel_param) plot(weeks, cdrsb_soc, type = "l", lwd = 2, xlab = "Weeks", ylab = "CDR-SB", main = "CDR-SB Progression for SOC", col = "grey", ylim = c(0,18), yaxt = "n") + lines(weeks, cdrsb_adu, lwd = 2, col = "red") + lines(weeks, cdrsb_leca, lwd = 2, col = "blue") + lines(weeks, cdrsb_dona, lwd = 2, col = "green") + axis(side = 2, at = seq(0,18,by =1)) ``` #2. Prepare empty matrices to recrod QALYs and Costs outcomes by cycle and by strategy: total_costs and total_QALYs have three empty data entries to match with the number of strategies (4) ```{r} # Sum costs and QALYs for each cycle at a time for each drug cycle_empty_array <- array(NA, dim = c(n_strategies, n_cycles), dimnames = list(treatment = t_names, cycle = 0:(n_cycles-1))) # per-cycle outcome template matrix cycle_state_costs <- cycle_trans_costs <- cycle_empty_array cycle_costs <- cycle_QALYs <- cycle_empty_array cycle_community_LY <- LE <- cycle_LYs <- cycle_empty_array # life expectancy; life-years cycle_QALE <- cycle_empty_array # quality-adjusted life expectancy cycle <- cycle_empty_array cycle_costs_dc <- cycle_LYs_dc <- cycle_QALYs_dc <- cycle_empty_array cycle_tx_cost <- cycle_costs_ontx <- cycle_LYs_ontx <- cycle_QALYs_ontx <- cycle_empty_array testtest<-cycle_empty_array testtest_origin<-cycle_empty_array #Total Cost & QALYs total_empty_array <- setNames(c(NA, NA, NA, NA), t_names) total_costs <- setNames(c(NA, NA, NA, NA), t_names) total_QALYs <- setNames(c(NA, NA, NA, NA), t_names) total_LYs <- setNames(c(NA, NA, NA, NA), t_names) total_ICER <- setNames(c(NA, NA, NA, NA), t_names) find_cycle <- setNames(c(NA, NA, NA, NA), t_names) total_tx_cost <- setNames(c(NA, NA, NA, NA), t_names) total_off_tx_tilModAD <- array(NA, dim = c(n_strategies, n_cycles), dimnames = list(treatment = t_names, cycle = 0:(n_cycles-1))) total_death <- total_off_tx <- total_on_tx <-total_off_tx_tilModAD #Sex-specific Total Cost & QALYs matrices total_costs_bymf <- array(NA, dim = c(n_strategies, n_sex), dimnames = list(Treatment = t_names, Sex = sex_2cat)) total_QALYs_bymf <- array(NA, dim = c(n_strategies, n_sex), dimnames = list(Treatment = t_names, Sex = sex_2cat)) #Stage-specific Total Cost & QALYs matrices total_costs_bystage <- array(NA, dim = c(n_strategies, n_states), dimnames = list(Treatment = t_names, Stage = s_names)) total_QALYs_bystage <- array(NA, dim = c(n_strategies, n_states), dimnames = list(Treatment = t_names, Stage = s_names)) #Setting-specific (community vs. ltc) Total Cost & QALYs matrices total_costs_bysetting <- array(NA, dim = c(n_strategies, n_setting), dimnames = list(Treatment = t_names, Setting = setting_names)) total_QALYs_bysetting <- array(NA, dim = c(n_strategies, n_setting), dimnames = list(Treatment = t_names, Setting = setting_names)) #Empty matrices for RESULTS df_ce <- data.frame(Strategy = t_names, Cost = numeric(n_strategies), QALY = numeric(n_strategies), LY = numeric(n_strategies), NMB = numeric(n_strategies), ICER = numeric(n_strategies), stringsAsFactors = FALSE) df_ce_2 <- df_ce ex_table<- data.frame(Strategy = t_names, "D/C before losing benefit" = numeric(n_strategies), "Continued tx until losing beneift" = numeric(n_strategies), "Any D/C" = numeric(n_strategies), "Total Death" = numeric(n_strategies), stringsAsFactors = FALSE) ###Breakdown analysis #cycle cycle_cost_pt_healthcare <- cycle_empty_array cycle_cost_cp_healthcare <- cycle_empty_array cycle_cost_pt_nonhealthcare <- cycle_empty_array cycle_cost_pt_oop <- cycle_empty_array cycle_cost_pt_productivity <- cycle_empty_array cycle_cost_cp_lostprod_institute <- cycle_empty_array cycle_cost_cp_lostprod_community <- cycle_empty_array cycle_ADE_cost <- cycle_empty_array cycle_mri_monitoring <- cycle_empty_array cycle_pt_ltc <- cycle_empty_array #total total_cost_pt_healthcare <- total_empty_array total_cost_cp_healthcare <- total_empty_array total_cost_pt_nonhealthcare <- total_empty_array total_cost_pt_oop <- total_empty_array total_cost_pt_productivity <- total_empty_array total_cost_cp_lostprod_institute <- total_empty_array total_cost_cp_lostprod_community <- total_empty_array total_ADE_cost <- total_empty_array total_mri_monitoring <- total_empty_array total_pt_ltc <- total_empty_array #Cost Check cycle_cost_check <- cycle_empty_array total_costs_check <- total_empty_array total_by_input <- data.frame(Strategy = t_names, total_tx_cost = numeric(n_strategies), total_cost_pt_healthcare = numeric(n_strategies), total_cost_cp_healthcare = numeric(n_strategies), total_cost_pt_nonhealthcare = numeric(n_strategies), total_cost_pt_oop = numeric(n_strategies), total_cost_pt_productivity = numeric(n_strategies), total_cost_cp_lostprod_institute = numeric(n_strategies), total_cost_cp_lostprod_community = numeric(n_strategies), total_ADE_cost = numeric(n_strategies), total_mri_monitoring = numeric(n_strategies), total_pt_ltc = numeric(n_strategies), stringsAsFactors = FALSE) ``` #3. Prepare population matrix that keeps track of health state transition by cycle ```{r} # Store population output for each cycle # state populations pop <- array(data = NA, dim = c(n_strategies, n_cycles), dimnames = list(treatment = t_names, cycle = 0:(n_cycles-1))) pop_check <- off_tx <- on_tx <- pop pop[,cycle = 1] <- n_cohort # pop["Mild AD", cycle = 1, ] <- 0 # pop["Moderate AD", cycle = 1, ] <- 0 # pop["Severe AD", cycle = 1, ] <- 0 # pop["Dead", cycle = 1, ] <- 0 # state population by sex pop_mf <- array(data = NA, dim = c(n_strategies, n_cycles, n_sex), dimnames = list(treatment = t_names, cycle = 0:(n_cycles-1), sex = sex_2cat)) pop_mf[ ,cycle = 1,"Male"] <- n_cohort*base_male_perc pop_mf[ ,cycle = 1,"Female"] <- n_cohort*(1-base_male_perc) #Empty Matrices to track On_Treatment and Off_Treatment population on_tx <- array(data = NA, dim = c(n_strategies, n_cycles), dimnames = list(treatment = t_names, cycle = 0:(n_cycles-1))) off_tx <- array(data = NA, dim = c(n_strategies,n_cycles), dimnames = list(treatment = t_names, cycle = 0:(n_cycles-1))) off_tx_bfModAD <- off_tx_tilModAD <- off_tx cycle_death <- array(data = NA, dim = c(n_strategies, n_cycles), dimnames = list(treatment = t_names, cycle = 0:(n_cycles-1))) cum_cycle_death <-cycle_death #Check # pop[cycle = 1, treatment = 4] # pop["Mild AD", cycle = 1, ] # pop["Moderate AD", cycle = 1, ] # pop["Severe AD", cycle = 1, ] # pop["Dead", cycle = 1, ] ``` #4. Simulate the cohort for the 4 strategies ```{r} #Run Simulation Model # for (k in 1:n_sex){ # sex <- if(k == 1) "Male" else if (k == 2) "Female" ARIA_tracker <- HA_tracker <- Fall_tracker <- IRR_tracker <-cum_cycle_death <- cycle_death <- modAD_cycle <- cycle_empty_array total_on_tx <- find_cycle <- setNames(c(NA, NA, NA, NA), t_names) total_off_tx_bfcomplete <- setNames(c(NA, NA, NA, NA), t_names) for (i in 1:n_strategies) {# outer loop: over strategies #Set initial age for this strategy age <- base_age male_perc <- base_male_perc female_perc <- base_female_perc for (j in 1:n_cycles){# inner loop: over cycles #Starting at cycle 2(1st week) ######################### Defining Discount Rates for Cost and Outcome ######################### dR_cost <- 1/(1+dc_costs)^(weeks[j]/weeksperyear) dR_outcome <- 1/(1+dc_outcomes)^(weeks[j]/weeksperyear) ######################### Calling CDR-SB & Using CDR-SB to determine stage ######################### cdrsb_options <- c( "1" = cdrsb_soc[j], "2" = cdrsb_adu[j], "3" = cdrsb_leca[j], "4" = cdrsb_dona[j] ) ######################### Define AD Stage based on CDR-SB ######################### stage <- ad_stage(cdrsb_options[i]) ######################### Calculate Population Survival in the current cycle ######################### if(j >= 2){ #Sex-specific population died in the current cycle cycle_mort_male <- mortality_combined_perwk(stage, age, "Male") * (pop[i,j-1]*male_perc) cycle_mort_female <- mortality_combined_perwk(stage, age, "Female") * (pop[i,j-1]*female_perc) #Population in the current cycle after removing mortality pop[i,j] <- pop[i,j-1] - sum(cycle_mort_male, cycle_mort_female, na.rm = TRUE) }else if(j == 1){ #Sex-specific population died in the current cycle cycle_mort_male <- 0 cycle_mort_female <- 0 pop[treatment = i, cycle = j] <- pop[treatment = i, cycle = j] } cycle_death[i, j] <- sum(cycle_mort_male, cycle_mort_female, na.rm = TRUE) cum_cycle_death[i,j] <-sum(cycle_death[i,1:j], na.rm=TRUE) ######################### TREATMENT DISCONTINUATION ######################### #1. D/c rate only applied for the first 6 months of starting the tx (ICER 2022 Beta-Amyloid CE model) dc_rate <- c( "1" = 0, "2" = dc_rate_6month_adu/weekspermonth, "3" = dc_rate_6month_leca/weekspermonth, "4" = dc_rate_6month_dona/weekspermonth ) #2. TRACKER FOR on-tx and off-tx pop w/ DC Rate #*Treatments are D/C when pts are in Moderate AD (lack of efficacy)*# #*D/C rate that is not due to AD stage is only applied in the first 6 months of starting tx*# if(i %in% c(2, 3, 4)){ if(j == 1){ off_tx[i, j] <- 0 #Everyone is on-treatment on_tx[i, j] <- pop[treatment = i, 1] }else if(j == 2){ off_tx[i,j] <- on_tx[i,1]*dc_rate[i] on_tx[i,j] <- pop[i,j]-off_tx[i,j] }else if(j >= 3 & j < weeksperyear/2){ off_tx[i,j] <- if(stage %in% c("MCI", "Mild AD")){ (off_tx[i,j-1] + (on_tx[i,j-1]*(dc_rate[i])))*(pop[i,j]/pop[i,j-1]) }else if(stage %in% c("Moderate AD", "Severe AD")){ off_tx[i, j-1]*((pop[i,j]/pop[i,j-1])) }else(stop("Error at d/c rate-1")) on_tx[i,j] <- pop[i,j]-off_tx[i,j] }else if(j >= 3 & j >= weeksperyear/2){ off_tx[i,j] <- off_tx[i,j-1]*(pop[i,j]/pop[i,j-1]) on_tx[i,j] <- pop[i,j]-off_tx[i,j] }else(stop("Error d/c rate-2")) }else if(i == 1){ on_tx[i, j] <- 0 off_tx[i, j] <- pop[i,j] }else(stop("Error d/c rate-3")) #2.1. D/C Treatment when reaching Moderate AD if(stage %in% c("Moderate AD", "Severe AD")){ on_tx[i,j] <- 0 off_tx[i,j] <- pop[i,j] } #3. Retrieving the CDR-SB trajectory that pts with treatment D/C follow cdrsb_dc <- kuhnel_cdrsb_dc(i,j)[j] stage_dc <- ad_stage(cdrsb_dc) ######################### ADVERSE EVENTS ######################### #*ASSUMED PT ONLY DEVELOPE ADE DURING FIRST 6 MONTHS OF THE TREATMENTS*# #*Define the week where ADE events stop*# ade_stop_week <- weeksperyear/2 #6 months = 24 weeks ##1. ARIA (Moderate, Symptomatic) (Pt resolves by 10th week and restart the therapy)## #Management: Monitoring MRI every 4 weeks for 12 weeks aria_tx <- c( "1" = 0, "2" = aria_pweek_adu, "3" = aria_pweek_leca, "4" = aria_pweek_dona ) ade_aria_pop <- if(j>1) on_tx[i,j-1]*aria_tx[i] else 0 ade_aria_cost <- if(j < ade_stop_week){#1 MD visit per event aria_mngm_cost*ade_aria_pop }else{0} ade_aria_disutility <- if(j < ade_stop_week){ aria_disutility*ade_aria_pop }else{0} #1.1 Place ARIA ppl to D/C until they are reintroduced after 12th week #aria_tx_dc <-aria_tx_restart if(j==1){ARIA_tracker[i,j]<-0} if(j < ade_stop_week & j > 1 & stage %in% c("MCI","Mild AD")){ off_tx[i,j] <- off_tx[i,j] + ade_aria_pop on_tx[i,j]<-pop[i,j]-off_tx[i,j] ARIA_tracker[i,j] <- ade_aria_pop }else if(j > (12+1)){ on_tx[i,j] <- on_tx[i,j]+sum(ARIA_tracker[i,j-12]*(pop[i,j]/pop[i,j-12]), na.rm=TRUE) off_tx[i,j] <- off_tx[i,j]-sum(ARIA_tracker[i,j-12]*(pop[i,j]/pop[i,j-12]), na.rm=TRUE) } #modAD_cycle <- cycle_empty_array modAD_cycle[i,j] <-if(stage == "MCI"){"MCI"}else if(stage == "Mild AD"){"Mild AD"}else if(stage == "Moderate AD"){"Moderate AD"}else if(stage == "Severe AD"){"Severe AD"} ##2. Headache## HA_tx <- c( "1" = 0, "2" = headache_pweek_adu, "3" = headache_pweek_leca, "4" = headache_pweek_dona ) ade_HA_pop <- if(j>1 & j < ade_stop_week){ on_tx[i,j-1]*HA_tx[i] }else{0} ade_HA_cost <- if(j < ade_stop_week){#1 MD visit per event drug_migraine_cost*ade_HA_pop }else{0} ade_HA_disutility <- if(j < ade_stop_week){ headache_disutility*ade_HA_pop }else{0} #disutility applied per event #HA_tracker <- cycle_empty_array HA_tracker[i,j] <- ade_HA_pop ##3. Fall (Mild-moderate ?)## fall_tx <- c( "1" = 0, "2" = fall_pweek_adu, "3" = fall_pweek_leca, "4" = fall_pweek_dona ) ade_fall_pop <- if(j >1 &j < ade_stop_week){ on_tx[i,j-1]*fall_tx[i] }else{0} ade_fall_cost <- if(j < ade_stop_week){#1 MD visit per event fall_mngm_cost*ade_fall_pop }else{0} ade_fall_disutility <- if(j < ade_stop_week){ fall_disutility*ade_fall_pop }else{0} #disutility applied per event #Fall_tracker <- cycle_empty_array Fall_tracker[i,j] <- ade_fall_pop ##4. IRR## irr_tx <- c( "1" = 0, "2" = irr_pweek_adu, "3" = irr_pweek_leca, "4" = irr_pweek_dona ) ade_irr_pop <- if(j > 1 & j < ade_stop_week){ on_tx[i,j-1]*irr_tx[i] }else{0} ade_irr_cost <- if(j < ade_stop_week){ irr_mngm_cost*ade_irr_pop }else{0} ade_irr_disutility <- if(j < ade_stop_week){ irr_disutility*ade_irr_pop }else{0} #disutility applied per event #IRR_tracker <- cycle_empty_array IRR_tracker[i,j] <- ade_irr_pop ADE_cost <- ade_aria_cost + ade_HA_cost + ade_fall_cost + ade_irr_cost ADE_disutilities <- ade_aria_disutility + ade_HA_disutility + ade_fall_disutility + ade_irr_disutility ######################### DRUG-RELATED COST ######################### if(i==1){ tx_pwk_cost <- 0 txadmin_pwk_cost <- 0 }else if(i == 2){ tx_pwk_cost <- if(j <= 9){ cost_txtitration_pmonth_adu/weekspermonth }else if(j <= 17){ (cost_txtitration_pmonth_adu*1.76)/weekspermonth }else if (j <= 25){ (cost_txtitration_pmonth_adu*3.53)/weekspermonth }else if (j > 25){ cost_tx_pmonth_adu/weekspermonth } txadmin_pwk_cost <- (cost_admin_pmonth_adu/weekspermonth) }else if(i == 3){ tx_pwk_cost <- (cost_tx_pmonth_leca/weekspermonth) txadmin_pwk_cost <- (cost_admin_pmonth_leca/weekspermonth) }else if(i == 4){ tx_pwk_cost <- if(j <= 13){ (cost_tx_pmonth_dona/2)/weekspermonth }else if(j > 13){ cost_tx_pmonth_dona/weekspermonth } txadmin_pwk_cost <- (cost_admin_pmonth_leca/weekspermonth) }else(stop("Error in drug cost")) ##Monitoring Cost for ROUTINE MRI for all pts in the beginning of their therapy #exclude pts developed ARIA as they will undergo separate MRI monitoring mri_monitoring <- monitor_adu(j, i)*(on_tx[i,j]-ade_aria_pop) ######################### COSTS ######################### #Off_Tx cycle_costs_dc[i, j] <- (off_tx[i,j]*((cost_pt_healthcare[stage_dc] + cost_cp_healthcare[stage_dc] + cost_pt_nonhealthcare[stage_dc] + cost_pt_oop[stage_dc] + cost_pt_productivity[stage_dc] + (ltc_percent[stage_dc]*(cost_cp_lostprod_institute[stage_dc] + cost_ltc[stage_dc])) + ((1-ltc_percent[stage_dc])*cost_cp_lostprod_community[stage_dc]))/weekspermonth))*dR_cost #On-Tx cycle_costs_ontx[i, j] <- ((on_tx[i,j]*(((cost_pt_healthcare[stage] + cost_cp_healthcare[stage] + cost_pt_nonhealthcare[stage] + cost_pt_oop[stage] + cost_pt_productivity[stage] + (ltc_percent[stage]*(cost_cp_lostprod_institute[stage] + cost_ltc[stage])) + ((1-ltc_percent[stage])*cost_cp_lostprod_community[stage]))/weekspermonth) +tx_pwk_cost+txadmin_pwk_cost))+ADE_cost+mri_monitoring)*dR_cost # # cycle_costs_ontx[i,j] <- (on_tx[i,j]*((cost_pt_healthcare[stage] + # cost_cp_healthcare[stage] + # cost_pt_nonhealthcare[stage] + # cost_pt_oop[stage] + # cost_pt_productivity[stage] + # (ltc_percent[stage]*(cost_cp_lostprod_institute[stage] + cost_ltc[stage])) + # ((1-ltc_percent[stage])*cost_cp_lostprod_community[stage]))/weekspermonth)+ # tx_pwk_cost + txadmin_pwk_cost + ADE_cost+mri_monitoring)*dR_cost #On_tx + off_tx Combined cycle_costs[i,j] <- cycle_costs_ontx[i,j] + cycle_costs_dc[i,j] ######################### SUMMING UTILITIES ######################### #off_tx cycle_LYs_dc[i,j] <- off_tx[i,j]/weeksperyear cycle_QALYs_dc[i,j]<-(cycle_LYs_dc[i, j]*((ltc_percent[stage_dc]*utility_ltc[stage_dc]) + ((1-ltc_percent[stage_dc])*utility_comm[stage_dc])-disutility_cp[stage_dc]))*dR_outcome #on_tx cycle_LYs_ontx[i,j] <- on_tx[i,j]/weeksperyear #cycle_community_LY_ontx[treatment = i, cycle = j] <- cycle_LYs[treatment = i, cycle = j]*(1-ltc_percent[stage]) cycle_QALYs_ontx[i,j] <- ((cycle_LYs_ontx[i,j] * ((ltc_percent[stage]*utility_ltc[stage]) + ((1-ltc_percent[stage])*utility_comm[stage]) - disutility_cp[stage]))- ADE_disutilities) * dR_outcome #Test testtest_origin[i,j]<-ADE_disutilities #On_tx + off_tx Combined cycle_LYs[treatment = i, cycle = j] <- cycle_LYs_ontx[treatment = i, cycle = j] + cycle_LYs_dc[treatment = i, cycle = j] cycle_QALYs[treatment = i, cycle = j] <- cycle_QALYs_ontx[treatment = i, cycle = j] + cycle_QALYs_dc[treatment = i, cycle = j] ######################### Updating Parameters for the next cycle ######################### #Update Age age <- base_age + (weeks[j]/weeksperyear) #Update Male and Female Percentage in the population if(j>=2){ male_perc <- ((pop[treatment = i, cycle = j-1]*male_perc)-cycle_mort_male)/pop[treatment = i, cycle = j] female_perc <- ((pop[treatment = i, cycle = j-1]*female_perc)-cycle_mort_female)/pop[treatment = i, cycle = j] }else{ male_perc <- male_perc female_perc <- female_perc } ######################### BreakDown Analysis ######################### # Patient Healthcare Costs cycle_cost_pt_healthcare[i,j] <- (((on_tx[i,j] * cost_pt_healthcare[stage]) + (off_tx[i,j] * cost_pt_healthcare[stage_dc])) / weekspermonth) * dR_cost # Caregiver Healthcare Costs cycle_cost_cp_healthcare[i,j] <- (((on_tx[i,j] * cost_cp_healthcare[stage]) + (off_tx[i,j] * cost_cp_healthcare[stage_dc])) / weekspermonth) * dR_cost # Patient Non-Healthcare Costs cycle_cost_pt_nonhealthcare[i,j] <- (((on_tx[i,j] * cost_pt_nonhealthcare[stage]) + (off_tx[i,j] * cost_pt_nonhealthcare[stage_dc])) / weekspermonth) * dR_cost # Patient Out-of-Pocket Costs cycle_cost_pt_oop[i,j] <- (((on_tx[i,j] * cost_pt_oop[stage]) + (off_tx[i,j] * cost_pt_oop[stage_dc])) / weekspermonth) * dR_cost # Patient Productivity Costs cycle_cost_pt_productivity[i,j] <- (((on_tx[i,j] * cost_pt_productivity[stage]) + (off_tx[i,j] * cost_pt_productivity[stage_dc])) / weekspermonth) * dR_cost # Caregiver Lost Productivity - Institutional Care cycle_cost_cp_lostprod_institute[i,j] <- (((on_tx[i,j] * cost_cp_lostprod_institute[stage] * ltc_percent[stage]) + (off_tx[i,j] * cost_cp_lostprod_institute[stage_dc] * ltc_percent[stage_dc])) / weekspermonth) * dR_cost # Caregiver Lost Productivity - Community Care cycle_cost_cp_lostprod_community[i,j] <- (((on_tx[i,j] * cost_cp_lostprod_community[stage] * (1 - ltc_percent[stage])) + (off_tx[i,j] * cost_cp_lostprod_community[stage_dc] * (1 - ltc_percent[stage_dc]))) / weekspermonth) * dR_cost # Long-Term Care Costs cycle_pt_ltc[i,j] <- (((on_tx[i,j] * cost_ltc[stage] * ltc_percent[stage]) + (off_tx[i,j] * cost_ltc[stage_dc] * ltc_percent[stage_dc])) / weekspermonth) * dR_cost # # Treatment Drug Cost (only for on-treatment patients) # cycle_cost_tx_drug[i,j] <- (on_tx[i,j] * tx_pwk_cost) * dR_cost # # # Treatment Administration Cost (only for on-treatment patients) # cycle_cost_tx_admin[i,j] <- (on_tx[i,j] * txadmin_pwk_cost) * dR_cost # Adverse Drug Event Costs (only for on-treatment patients) cycle_ADE_cost[i,j] <- ADE_cost* dR_cost # MRI Monitoring Costs (only for on-treatment patients) cycle_mri_monitoring[i,j] <- mri_monitoring * dR_cost #Cycle Tx Cost cycle_tx_cost[i,j] <- (on_tx[i,j]*(tx_pwk_cost+txadmin_pwk_cost))*dR_cost # #Treatment Cost: drug cost + administration cost # cycle_tx_cost[i,j] <- (on_tx[i,j]*(tx_pwk_cost+txadmin_pwk_cost))*dR_cost # #Check cycle total cost cycle_cost_check[i,j] <- cycle_tx_cost[i,j] + cycle_cost_pt_healthcare[i,j] + cycle_cost_cp_healthcare[i,j] + cycle_cost_pt_nonhealthcare[i,j] + cycle_cost_pt_oop[i,j] + cycle_cost_pt_productivity[i,j] + cycle_cost_cp_lostprod_institute[i,j] + cycle_cost_cp_lostprod_community[i,j] + cycle_pt_ltc[i,j] + cycle_ADE_cost[i,j] + cycle_mri_monitoring[i,j] #Population Check pop_check[i,j] <- on_tx[i,j] + off_tx[i,j] + cum_cycle_death[i, j] } #Total Cost & QALYs total_costs[i] <- sum(cycle_costs[treatment = i, ],na.rm=TRUE) total_QALYs[i] <- sum(cycle_QALYs[treatment = i, ],na.rm=TRUE) total_LYs[i] <- sum(cycle_LYs[treatment = i, ], na.rm=TRUE) total_ICER[i] <- (total_costs[i]-total_costs[1])/(total_QALYs[i]-total_QALYs[1]) df_ce[i, c("Cost")] <- as.numeric(total_costs[i]) df_ce[i, c("QALY")] <- as.numeric(total_QALYs[i]) df_ce[i, c("LY")] <- as.numeric(total_LYs[i]) df_ce[i, c("NMB")] <- (as.numeric(total_QALYs[i])*200000) - as.numeric(total_costs[i]) df_ce[i, c("ICER")] <- as.numeric(total_ICER[i]) #ADDITIONAL ANALYSIS RESULTS find_cycle[i] <- cycle[treatment=i,1] total_off_tx_bfcomplete[i] <- as.numeric(off_tx[i, find_cycle[i]-1]) # as.numeric(off_tx_bfModAD[treatment = i, 1]) total_on_tx[i] <- as.numeric(on_tx[treatment=i, find_cycle[i]-1]) total_off_tx[i] <- off_tx[treatment=i, n_cycles] total_death[i] <- sum(cycle_death[treatment = i,], na.rm=TRUE) ex_table[i,c(2,3,4,5)] <- c(total_off_tx_bfcomplete[i], total_on_tx[i], total_off_tx[i], total_death[i]) #Treatment Cost (Drug Cost + Administration Cost) total_tx_cost[i] <- sum(cycle_tx_cost[i, ], na.rm=TRUE) ###BreakDown Analysis #total total_cost_pt_healthcare[i] <- sum(cycle_cost_pt_healthcare[treatment = i, ],na.rm=TRUE) total_cost_cp_healthcare[i] <- sum(cycle_cost_cp_healthcare[treatment = i, ],na.rm=TRUE) total_cost_pt_nonhealthcare[i] <- sum(cycle_cost_pt_nonhealthcare[treatment = i, ],na.rm=TRUE) total_cost_pt_oop[i] <- sum(cycle_cost_pt_oop[treatment = i, ],na.rm=TRUE) total_cost_pt_productivity[i] <- sum(cycle_cost_pt_productivity[treatment = i, ],na.rm=TRUE) total_cost_cp_lostprod_institute[i] <- sum(cycle_cost_cp_lostprod_institute[treatment = i, ],na.rm=TRUE) total_cost_cp_lostprod_community[i] <- sum(cycle_cost_cp_lostprod_community[treatment = i, ],na.rm=TRUE) total_ADE_cost[i] <- sum(cycle_ADE_cost[treatment = i, ],na.rm=TRUE) total_mri_monitoring[i] <- sum(cycle_mri_monitoring[treatment = i, ],na.rm=TRUE) total_pt_ltc[i] <- sum(cycle_pt_ltc[treatment = i, ],na.rm=TRUE) total_by_input[i, c("total_tx_cost")] <- as.numeric(total_tx_cost[i]) total_by_input[i, c("total_cost_pt_healthcare")] <- as.numeric(total_cost_pt_healthcare[i]) total_by_input[i, c("total_cost_cp_healthcare")] <- as.numeric(total_cost_cp_healthcare[i]) total_by_input[i, c("total_cost_pt_nonhealthcare")] <- as.numeric(total_cost_pt_nonhealthcare[i]) total_by_input[i, c("total_cost_pt_oop")] <- as.numeric(total_cost_pt_oop[i]) total_by_input[i, c("total_cost_pt_productivity")] <- as.numeric(total_cost_pt_productivity[i]) total_by_input[i, c("total_cost_cp_lostprod_institute")] <- as.numeric(total_cost_cp_lostprod_institute[i]) total_by_input[i, c("total_cost_cp_lostprod_community")] <- as.numeric(total_cost_cp_lostprod_community[i]) total_by_input[i, c("total_ADE_cost")] <- as.numeric(total_ADE_cost[i]) total_by_input[i, c("total_mri_monitoring")] <- as.numeric(total_mri_monitoring[i]) total_by_input[i, c("total_pt_ltc")] <- (as.numeric(total_pt_ltc[i])) #Total Cost Check total_costs_check[i]<-sum(cycle_cost_check[i,], na.rm=TRUE) } ``` ########After running the simulation loop above, cross-check to make sure it was performed correctly.######### ```{r} #Check Population print(pop_check[4,]) # # # total_tx_cost <- total_costs # cycle_tx_cost <- cycle_costs # cycle<- array(data = NA, # dim = c(n_strategies, n_cycles), # dimnames = list(treatment = t_names, # cycle = 0:(n_cycles-1))) # #Gathering data for Supplemental Material # #A. Track AD Stage Transitions # write.csv(modAD_cycle, "modAD_cycle_2024_raw.csv", row.names = FALSE) # #B. Track people died before each stage (refer to A. for transition time) # write.csv(cycle_death, "cycle_death_2024_raw.csv", row.names = FALSE) # #C. Track people D/C'ed from Tx # write.csv(off_tx, "Off_Tx_2024_raw.csv", row.names = FALSE) # #D. Track Adverse event incidence rate # write.csv(ARIA_tracker, "ARIA_tracker_2024_raw.csv", row.names = FALSE) # write.csv(HA_tracker, "HA_tracker_2024_raw.csv", row.names = FALSE) # write.csv(Fall_tracker, "Fall_tracker_2024_raw.csv", row.names = FALSE) # write.csv(IRR_tracker, "IRR_tracker_2024_raw.csv", row.names = FALSE) # write.csv(pop, "Pop_raw.csv") sum(cycle_costs_ontx[2,],na.rm=TRUE) sum(cycle_costs_dc[2,],na.rm=TRUE) print(ex_table) print(df_ce) # print(ex_table) print(total_by_input) write.csv(total_by_input, "Total_by_input.csv", row.names = FALSE) #Check Costs sum(cycle_costs[1, ],na.rm=TRUE) sum(total_by_input[1,c(2:12)],na.rm=TRUE) # write.csv(pars_basecase, "pars_basecase.csv", row.names = FALSE) ``` ############################ RESULTS ############################ #Table 2A: Comparative Cost-Effectiveness Results 9. Check the total cost and qaly by strategy. Calculate incremental cost and qaly and ICER *Four strategies are in the increasing order in both cost and qaly. Hence, no strategy is strongly dominated. (LEAGUE TABLE)* ```{r} #Make League Table write.csv(df_ce, "CE_4strategies_raw.csv", row.names = FALSE) df_ce icer_mab <- calculate_icers(cost = df_ce[,c("Cost")], effect = df_ce[,c("QALY")], strategies = t_names) icer_mab <- icer_mab %>% arrange(Cost) icer_mab icer_mab %>% kable() %>% kable_styling() plot(icer_mab, label = "all") + theme_classic() + ggtitle("Cost-effectiveness of Monoclonal Antibody Treatments for AD") ``` `