Recalibrated results and code cleaning

SA_saturation.R

@@ -10,44 +10,54 @@ rm(list=ls())
 args = commandArgs(trailingOnly=TRUE)
 
 ## Model setup parameters ##
-yr.first     <- 2016        # starting year of simulation 
-yr.last      <- 2022        # end year of simulation (also the year for evaluation)
-d.c          <- 0.03        # discounting of costs by 3%
-seed         <- 2021
-sw.EMS.ODloc <- "ov"  #Please choose from "ov" (using average overall) or "sp" (region-specific) for overdose setting parameter, default is "ov"
+yr_start <- 2016        # starting year of simulation 
+yr_end <- 2022        # end year of simulation (also the year for evaluation)
+discount.rate <- 0.03        # discounting of costs by 3%
+# seed <- 2021
 
 ## LOAD packages and functions
+library("argparser")
 library(openxlsx)
 library(dplyr)
 library(abind)
 library(tictoc)
-source("Profound-Function-PopInitialization.R")
-source("Profound-Function-TransitionProbability.R")
-source("Profound-Function-Microsimulation.R")
-source("Profound-DecisionTree.R")
-source("Profound-DataInput.R")
-source("Profound-Function-NxAvailAlgm.R")
-source("Profound-CEA.R")
-source("Profound-InputOutput-Setup.R")
+source("transition_probability.R")
+source("microsim.R")
+source("decision_tree.R")
+source("data_input.R")
+source("naloxone_availability.R")
+source("cost_effectiveness.R")
 
-## Initialize the study population - people who are at risk of opioid overdose
-pop.info  <- c("sex", "race", "age", "residence", "curr.state", "OU.state", "init.age", "init.state", "ever.od", "fx")
-if(file.exists(paste0("Inputs/InitialPopulation.rds"))){
-  init.pop  <- readRDS(paste0("Inputs/InitialPopulation.rds"))
-} else if (!file.exists(paste0("Inputs/InitialPopulation.rds"))){
-  init.pop  <- pop.initiation(initials = initials, seed=seed)
-  saveRDS(init.pop, paste0("Inputs/InitialPopulation.rds"))
-}
+args <- arg_parser("arguments")
+args <- add_argument(args, "--seed", help = "seed for random numbers", default = 2021)
+args <- add_argument(args, "--regional", help = "flag to run regional model", flag = TRUE)
+args <- add_argument(args, "--outfile", help = "file to store outputs", default = "OverdoseDeath_RIV1_0.csv")
+args <- add_argument(args, "--ppl", help = "file with initial ppl info", default = "Inputs/init_pop.rds")
+argv <- parse_args(args)
+seed <- as.integer(argv$seed)
+init_ppl.file <- argv$ppl
+source("io_setup.R")
+
+# ## Initialize the study population - people who are at risk of opioid overdose
+# pop.info  <- c("sex", "race", "age", "residence", "curr.state", "OU.state", "init.age", "init.state", "ever.od", "fx")
+# if(file.exists(paste0("Inputs/InitialPopulation.rds"))){
+#   init.pop  <- readRDS(paste0("Inputs/InitialPopulation.rds"))
+# } else if (!file.exists(paste0("Inputs/InitialPopulation.rds"))){
+#   init.pop  <- pop.initiation(initials = initials, seed=seed)
+#   saveRDS(init.pop, paste0("Inputs/InitialPopulation.rds"))
+# }
 
 ## Load calibrated parameters and seeds
 sim.data.ls <- readRDS(file = paste0("calibration/CalibratedData.rds"))
 sim.seed    <- readRDS(file = paste0("calibration/CalibratedSeed.rds"))
-sim.seed    <- sim.seed[1:50]
+sim.seed    <- sim.seed[1:100]
 
 #define different statewide program expansion scenarios, including a 0 level and a saturation level
-scenario.name <- c("Zero", "Status Quo", "10,000", "20,000", "50,000", "100,000", "150,000", "200,000", "250,000", "300,000", "350,000", "400,000", "450,000", "500,000", "600,000", "700,000", "800,000", "900,000", "1M")
+# scenario.name <- c("Zero", "Status Quo", "10,000", "20,000", "30,000", "40,000", "50,000", "60,000", "70,000", "80,000", "90,000", "100,000")
 sq.nlx     <- 8241
-expand.nlx <- c(0, 8241, 10000, 20000, 50000, 100000, 150000, 200000, 250000, 300000, 350000, 400000, 450000, 500000, 600000, 700000, 800000, 900000, 1000000)
+# expand.nlx <- c(0, 8241, 10000, 20000, 30000, 40000, 50000, 60000, 70000, 80000, 90000, 100000)
+scenario.name <- c("Zero", "Status Quo", "200,000", "400,000","600,000", "800,000", "1,000,000")
+expand.nlx <- c(0, 8241, 200000, 400000, 600000, 800000, 1000000)
 expand.level  <- expand.nlx / sq.nlx
 od.death.mx.last <- od.death.mx.totl <- od.death.mx.wtns <- od.death.mx.wtns.tl <- matrix(0, nrow = length(sim.seed), ncol = length(scenario.name))
 colnames(od.death.mx.last) <- colnames(od.death.mx.totl) <- colnames(od.death.mx.wtns) <- colnames(od.death.mx.wtns.tl) <- scenario.name
@@ -56,33 +66,33 @@ for (ss in 1:length(sim.seed)){
   # tic()
   print(paste0("Parameter set: ", ss))
   vparameters.temp<- sim.data.ls[[ss]]
-  vparameters.temp$NxDataPharm$pe  <- 0 # ATTN: This is temporary, used to manually remove pharamacy naloxone
-  vparameters.temp$mor_Nx <- vparameters.temp$mor_bl * (1-0.95) # ATTN: This is temporary, used to manually force naloxone effect to be 90% (reducing probability of dying from an overdose)
-  sim_sq          <- MicroSim(init.pop, vparameters = vparameters.temp, n.t, v.state, d.c, PT.out = FALSE, Str = "SQ", seed = sim.seed[ss])        # run for status quo
-  od.death.mx.last[ss, "Status Quo"]    <- sum(sim_sq$m.oddeath[(n.t-11):n.t, ])
-  od.death.mx.totl[ss, "Status Quo"]    <- sum(sim_sq$m.oddeath[(n.t-35):n.t, ])
-  od.death.mx.wtns[ss, "Status Quo"]    <- sum(sim_sq$v.oddeath.w[(n.t-11):n.t])
-  od.death.mx.wtns.tl[ss, "Status Quo"] <- sum(sim_sq$v.oddeath.w[(n.t-35):n.t])
+  # vparameters.temp$NxDataPharm$pe  <- 0 # ATTN: This is temporary, used to manually remove pharamacy naloxone
+  # vparameters.temp$mor_Nx <- vparameters.temp$mor_bl * (1-0.95) # ATTN: This is temporary, used to manually force naloxone effect to be 90% (reducing probability of dying from an overdose)
+  sim_sq <- MicroSim(init_ppl, params = vparameters.temp, timesteps, agent_states, discount.rate, PT.out = FALSE, strategy = "SQ", seed = sim.seed[ss])        # run for status quo
+  od.death.mx.last[ss, "Status Quo"]    <- sum(sim_sq$m.oddeath[(timesteps-11):timesteps, ])
+  od.death.mx.totl[ss, "Status Quo"]    <- sum(sim_sq$m.oddeath[(timesteps-35):timesteps, ])
+  od.death.mx.wtns[ss, "Status Quo"]    <- sum(sim_sq$v.oddeath.w[(timesteps-11):timesteps])
+  od.death.mx.wtns.tl[ss, "Status Quo"] <- sum(sim_sq$v.oddeath.w[(timesteps-35):timesteps])
 
   exp.lv <- 0
-  sim_pg  <- MicroSim(init.pop, vparameters = vparameters.temp, n.t, v.state, d.c, PT.out = FALSE, Str = "expand", seed = sim.seed[ss]) # run for program scenario
-  od.death.mx.last[ss, "Zero"]    <- sum(sim_pg$m.oddeath[(n.t-11):n.t, ])
-  od.death.mx.totl[ss, "Zero"]    <- sum(sim_pg$m.oddeath[(n.t-35):n.t, ])
-  od.death.mx.wtns[ss, "Zero"]    <- sum(sim_pg$v.oddeath.w[(n.t-11):n.t])
-  od.death.mx.wtns.tl[ss, "Zero"] <- sum(sim_pg$v.oddeath.w[(n.t-35):n.t])
+  sim_pg <- MicroSim(init_ppl, params = vparameters.temp, timesteps, agent_states, discount.rate, PT.out = FALSE, strategy = "expand", seed = sim.seed[ss]) # run for program scenario
+  od.death.mx.last[ss, "Zero"]    <- sum(sim_pg$m.oddeath[(timesteps-11):timesteps, ])
+  od.death.mx.totl[ss, "Zero"]    <- sum(sim_pg$m.oddeath[(timesteps-35):timesteps, ])
+  od.death.mx.wtns[ss, "Zero"]    <- sum(sim_pg$v.oddeath.w[(timesteps-11):timesteps])
+  od.death.mx.wtns.tl[ss, "Zero"] <- sum(sim_pg$v.oddeath.w[(timesteps-35):timesteps])
 
   for (jj in 3:length(expand.level)){
     exp.lv  <- expand.level[jj]
-    sim_pg  <- MicroSim(init.pop, vparameters = vparameters.temp, n.t, v.state, d.c, PT.out = FALSE, Str = "expand", seed = sim.seed[ss]) # run for program scenario
-    od.death.mx.last[ss, jj]    <- sum(sim_pg$m.oddeath[(n.t-11):n.t, ])
-    od.death.mx.totl[ss, jj]    <- sum(sim_pg$m.oddeath[(n.t-35):n.t, ])
-    od.death.mx.wtns[ss, jj]    <- sum(sim_pg$v.oddeath.w[(n.t-11):n.t])
-    od.death.mx.wtns.tl[ss, jj] <- sum(sim_pg$v.oddeath.w[(n.t-35):n.t])
+    sim_pg <- MicroSim(init_ppl, params = vparameters.temp, timesteps, agent_states, discount.rate, PT.out = FALSE, strategy = "expand", seed = sim.seed[ss]) # run for program scenario
+    od.death.mx.last[ss, jj]    <- sum(sim_pg$m.oddeath[(timesteps-11):timesteps, ])
+    od.death.mx.totl[ss, jj]    <- sum(sim_pg$m.oddeath[(timesteps-35):timesteps, ])
+    od.death.mx.wtns[ss, jj]    <- sum(sim_pg$v.oddeath.w[(timesteps-11):timesteps])
+    od.death.mx.wtns.tl[ss, jj] <- sum(sim_pg$v.oddeath.w[(timesteps-35):timesteps])
   }
   # toc()
 }
 
-write.xlsx(od.death.mx.last, file = ("Ignore/Saturation/Saturation_3Y_TotalODdeaths_last.xlsx"), row.names = F)
-write.xlsx(od.death.mx.totl, file = ("Ignore/Saturation/Saturation_3Y_TotalODdeaths_total.xlsx"), row.names = F)
-write.xlsx(od.death.mx.wtns, file = ("Ignore/Saturation/Saturation_3Y_WitnessedDeaths_last.xlsx"), row.names = F)
-write.xlsx(od.death.mx.wtns.tl, file = ("Ignore/Saturation/Saturation_3Y_WitnessedDeaths_total.xlsx"), row.names = F)
+write.xlsx(od.death.mx.last, file = ("Ignore/Saturation/Saturation_3Y_TotalODdeaths_last2.xlsx"), row.names = F)
+write.xlsx(od.death.mx.totl, file = ("Ignore/Saturation/Saturation_3Y_TotalODdeaths_total2.xlsx"), row.names = F)
+write.xlsx(od.death.mx.wtns, file = ("Ignore/Saturation/Saturation_3Y_WitnessedDeaths_last2.xlsx"), row.names = F)
+write.xlsx(od.death.mx.wtns.tl, file = ("Ignore/Saturation/Saturation_3Y_WitnessedDeaths_total2.xlsx"), row.names = F)

analyze_calibration.R

@@ -30,15 +30,33 @@ source("data_input.R")
 
 calibration_results <- NULL
 
-for (batch.ind in 1:5) {
+for (batch.ind in 1:56) {
+  temp_results <- readRDS(paste0("calibration/CalibrationOutputs", batch.ind, ".rds"))
+  calibration_results <- rbind(calibration_results, temp_results)
+}
+temp_results[,] <- 0
+calibration_results <- rbind(calibration_results, temp_results)
+for (batch.ind in 58:70) {
+  temp_results <- readRDS(paste0("calibration/CalibrationOutputs", batch.ind, ".rds"))
+  calibration_results <- rbind(calibration_results, temp_results)
+}
+temp_results[,] <- 0
+calibration_results <- rbind(calibration_results, temp_results)
+for (batch.ind in 72:75) {
+  temp_results <- readRDS(paste0("calibration/CalibrationOutputs", batch.ind, ".rds"))
+  calibration_results <- rbind(calibration_results, temp_results)
+}
+temp_results[,] <- 0
+calibration_results <- rbind(calibration_results, temp_results)
+for (batch.ind in 77:100) {
   temp_results <- readRDS(paste0("calibration/CalibrationOutputs", batch.ind, ".rds"))
   calibration_results <- rbind(calibration_results, temp_results)
 }
 rm(temp_results)
 
 calibration_params <- readRDS(paste0("calibration/prep_calibration_data/Calib_par_table.rds"))
 
-calibration_results <- cbind(calibration_results, calibration_params)
+calibration_results <- cbind(calibration_results, calibration_params[1:dim(calibration_results)[1], ])
 
 # read in workbook
 WB <- loadWorkbook("Inputs/MasterTable.xlsx")
@@ -65,11 +83,11 @@ for (ss in 1:nrow(calibration_results)) {
   calibration_results[ss, "gof"] <- gof
 }
 
-# sort by goodness of fit and select top 100 fits
+# sort by goodness of fit and select top 1000 fits
 sorted.mx <- calibration_results[order(calibration_results[, "gof"], decreasing = F), ]
-cal_sample <- 40
+cal_sample <- 1000
 calibration_results_subset <- sorted.mx[1:cal_sample, ]
-write.xlsx(calibration_results,
+write.xlsx(calibration_results_subset,
   file = paste0("calibration/Calibrated_results.xlsx"),
   col.names = T, row.names = F
 )
@@ -122,6 +140,14 @@ for (cc in 1:nrow(calibration_results_subset)) {
 saveRDS(calibrated_parameters, file = paste0("calibration/CalibratedData.rds"))
 saveRDS(calibration_results_subset[, "seed"], file = paste0("calibration/CalibratedSeed.rds"))
 
+#### plot calibrated results
+cal_sample <- 500
+calibration_results_subset <- read.xlsx( "calibration/Calibrated_results.xlsx")
+calibration_results_subset <- calibration_results_subset[1:cal_sample, ]
+# read in workbook
+WB <- loadWorkbook("Inputs/MasterTable.xlsx")
+Target <- read.xlsx(WB, sheet = "Target")
+tar.data <- Target$pe
 
 
 ## Plot ##
@@ -162,7 +188,7 @@ abline(h = 0)
 mean_fxdeath <- apply(calibration_results_subset[, c("fx.death16", "fx.death17", "fx.death18", "fx.death19")], 2, median)
 # REVIEWED: change from hardcoded ylim; mean instead of median for everything
 plot(
-  x = 2016:2019, tar.data[5:8], col = "black", pch = 18, xlab = "Year", ylab = "Proportion of OOD deaths with fentanyl present", cex = 1.2, cex.axis = 1.2, cex.lab = 1.3,
+  x = 2016:2019, tar.data[5:8], col = "black", pch = 18, xlab = "Year", ylab = "Percentage of opioid overdose deaths involving fentanyl", cex = 1.2, cex.axis = 1.2, cex.lab = 1.3,
   xaxt = "n", yaxt = "n", ylim = c(0, 1), frame.plot = FALSE
 )
 # title("A", adj = 0, line =-0.5)
@@ -192,7 +218,7 @@ axis(1, at = 2016:2019, pos = 0, lwd.ticks = 0, cex.axis = 1.2)
 abline(h = 0)
 
 par(mfrow = c(1, 1))
-legend("bottom",
+legend("top",
   legend = c("Target", "Model: mean", "Model: simulation"),
   col = c("black", "red", adjustcolor("indianred1", alpha = 0.2)),
   pch = c(16, NA, NA),