diff --git a/Inputs/MasterTable.xlsx b/Inputs/MasterTable.xlsx index 818673b..1469ba0 100644 Binary files a/Inputs/MasterTable.xlsx and b/Inputs/MasterTable.xlsx differ diff --git a/Inputs/ProgramData.xlsx b/Inputs/ProgramData.xlsx new file mode 100644 index 0000000..e050c62 Binary files /dev/null and b/Inputs/ProgramData.xlsx differ diff --git a/SA_program_time.R b/SA_program_time.R deleted file mode 100644 index 141fbbd..0000000 --- a/SA_program_time.R +++ /dev/null @@ -1,95 +0,0 @@ -####################################################################################################### -##################### Sensitivity analysis - statewide program ################# -####################################################################################################### -# To assess the impact of statewide program (naloxone increase by the same multiplier in all regions) - -############################################################################# -# 1. SET directory and workspace -############################################################################# -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" - -# Program data -library(openxlsx) -pg.data <- read.xlsx("Ignore/OEND_program.xlsx", sheet = "Project Weber") #load data for the program, including high/low-risk and last year's numbers -pg.levels <- c(2, 5, 10, 20, 50) #define program expansion levels -scenario.name <- c("Status Quo", "Double", "5 times", "10 times", "20 times", "50 times") -pg.add.array <- array(0, dim = c(length(pg.levels), 2, dim(pg.data)[1])) # array for the added naloxone kits (only acount for the additional kits) due to program expansion, stratified by high/low-risk program and regions -dimnames(pg.add.array)[[2]] <- c("high", "low") -for (i in 1:length(pg.levels)){ - if (pg.data$Risk[1] == "high"){ - pg.add.array[i, "high", ] <- round(pg.data$Volume[1] * pg.data$Proportion * (pg.levels[i]-1), 0) - pg.add.array[i, "low", ] <- 0 - } else { - pg.add.array[i, "high", ] <- 0 - pg.add.array[i, "low", ] <- round(pg.data$Volume[1] * pg.data$Proportion * (pg.levels[i]-1), 0) - } -} - -## LOAD packages and functions -library(openxlsx) -library(dplyr) -library(abind) -source("Profound-Function-PopInitialization.R") -source("Profound-Function-TransitionProbability.R") -source("Profound-Function-Microsimulation(Trial).R") -source("Profound-DecisionTree.R") -source("Profound-DataInput.R") -source("Profound-Function-NxAvailAlgm.R") -source("Profound-CEA.R") -source("Profound-InputOutput-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] - -m.oddeath.preb <- m.oddeath.il.lr <- m.oddeath.il.hr <- m.nlx.mn <- rep(0, times = n.t) -od.death.mx.first <- od.death.mx.total <- od.death.mx.last <- matrix(0, nrow = length(sim.seed), ncol = 1+length(pg.levels)) -colnames(od.death.mx.first) <- colnames(od.death.mx.total) <- colnames(od.death.mx.last) <- scenario.name -v.nlx.mn <- matrix(0, nrow= n.t, ncol = length(scenario.name)) -colnames(v.nlx.mn) <- scenario.name -for (ss in 1:length(sim.seed)){ - 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.9) # 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 - v.nlx.mn[ , "Status Quo"] <- sim_sq$m.nlx.mn - od.death.mx.first[ss, "Status Quo"] <- sum(sim_sq$m.oddeath[49:60, ]) - od.death.mx.last[ss, "Status Quo"] <- sum(sim_sq$m.oddeath[(n.t-11):n.t, ]) - od.death.mx.total[ss, "Status Quo"] <- sum(sim_sq$m.oddeath[49:n.t, ]) - - for (ll in 1:dim(pg.add.array)[1]){ - vparameters.temp$pg.add <- pg.add.array[ll, , ] - sim_pg <- MicroSim(init.pop, vparameters = vparameters.temp, n.t, v.state, d.c, PT.out = FALSE, Str = "program", seed = sim.seed[ss]) # run for program scenario - v.nlx.mn[ , ll+1] <- sim_pg$m.nlx.mn - od.death.mx.first[ss, ll+1] <- sum(sim_pg$m.oddeath[49:60, ]) - od.death.mx.last[ss, ll+1] <- sum(sim_pg$m.oddeath[(n.t-11):n.t, ]) - od.death.mx.total[ss, ll+1] <- sum(sim_pg$m.oddeath[49:n.t, ]) - } -} -write.xlsx(od.death.mx.first, file="Ignore/SA/time/Program3Yfirst.xlsx", - col.names = T, row.names = F) -write.xlsx(od.death.mx.total, file="Ignore/SA/time/Program3Ytotal.xlsx", - col.names = T, row.names = F) -write.xlsx(od.death.mx.last, file="Ignore/SA/time/Program3Ylast.xlsx", - col.names = T, row.names = F) -write.xlsx(v.nlx.mn, file="Ignore/SA/time/NaloxoneMonthly.xlsx", - col.names = T, row.names = F) diff --git a/SA_saturation.R b/SA_saturation.R index f31e4b4..a6196b2 100644 --- a/SA_saturation.R +++ b/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) \ No newline at end of file +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) \ No newline at end of file diff --git a/analyze_calibration.R b/analyze_calibration.R index 83ef436..427df72 100644 --- a/analyze_calibration.R +++ b/analyze_calibration.R @@ -30,7 +30,25 @@ 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) } @@ -38,7 +56,7 @@ 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), diff --git a/calibration/CalibratedData.rds b/calibration/CalibratedData.rds index 20c0c42..f154040 100644 Binary files a/calibration/CalibratedData.rds and b/calibration/CalibratedData.rds differ diff --git a/calibration/CalibratedSeed.rds b/calibration/CalibratedSeed.rds index 32fc9ad..d99b111 100644 Binary files a/calibration/CalibratedSeed.rds and b/calibration/CalibratedSeed.rds differ diff --git a/calibration/Calibrated_results.xlsx b/calibration/Calibrated_results.xlsx index 0fef4e7..54b8a9d 100644 Binary files a/calibration/Calibrated_results.xlsx and b/calibration/Calibrated_results.xlsx differ diff --git a/main.R b/main.R index 19f40fa..1883842 100755 --- a/main.R +++ b/main.R @@ -1,21 +1,20 @@ #!/usr/bin/env Rscript ############################################################################################### -###################### PROFOUND Naloxone Distribution model #### 2020 ######################### +###################### PROFOUND Naloxone Distribution model #### 2022 ######################### ############################################################################################### # Main module for the microsimulation of the Profound Naloxone distribution model: # # Author: Xiao Zang, PhD; Shayla Nolen, MPH; Sam Bessy, MSc # Marshall Lab, Department of Epidemiology, Brown University # Created: May 06, 2020 -# Last update: July 17, 2021 # ATTN: currently not used but will incorporate other main analysis modules here later ############################################################################################### ############################################################################################### #### Individual-based microsimulation #### -#### 7 health states: prescribed (preb), unregulated-injection (unreg.inj) #### -#### unregulated-noninjection (unreg.nin) #### +#### 7 health states: prescribed (preb), illicit-injection (il.hr) - high-risk #### +#### illicit-noninjection (il.lr) - low-risk #### #### inactive (inact), non-opioid drug use (NODU) - stimulant, #### #### relapsed (relap), death (dead) #### #### 1 health event: Overdose #### @@ -30,7 +29,7 @@ # 1. SET directory and workspace ############################################################################# rm(list = ls()) -library("argparser") +library(argparser) library(dplyr) library(tictoc) library(openxlsx) @@ -61,10 +60,13 @@ source("io_setup.R") ############## RI modeling analysis: distributing 10,000 kits through different programs ################## ########################################################################################################### # Define different program scenarios for distributing additional 10,000 kits and initialize matrices and arrrays for results -scenario.name <- c("Status Quo", "SSP_10K", "StreetOutreach_10K", "MailEvent_10K", "Healthcare_10K") -mx.od.death.last <- mx.od.death.totl <- mx.od.death.wtns.last <- mx.od.death.wtns.totl <- matrix(0, nrow = length(sim.seed), ncol = length(scenario.name)) +sim.seed <- sim.seed[1:500] +scenario.name <- c("Status Quo", "SSP_10K", "MailEvent_10K", "Healthcare_10K") +mx.od.death.last <- matrix(0, nrow = length(sim.seed), ncol = length(scenario.name)) mx.costs.totl <- matrix(0, nrow = length(sim.seed), ncol = length(scenario.name)) -colnames(mx.od.death.last) <- colnames(mx.od.death.totl) <- colnames(mx.od.death.wtns.last) <- colnames(mx.od.death.wtns.totl) <- colnames(mx.costs.totl) <-scenario.name +colnames(mx.od.death.last) <- colnames(mx.costs.totl) <-scenario.name +mx.od.death.1st <- mx.od.death.2nd <- mx.od.death.totl <- mx.od.death.wtns.1st <- mx.od.death.wtns.2nd <- mx.od.death.wtns.last <- mx.od.death.wtns.totl <- mx.od.death.last + array.od.death.rgn <- array(0, dim = c(num_regions, length(scenario.name), length(sim.seed))) array.oend.nlx.rgn <- array(0, dim = c(num_regions, length(scenario.name), length(sim.seed))) dimnames(array.od.death.rgn)[[1]] <- dimnames(array.oend.nlx.rgn)[[1]] <- v.region @@ -74,77 +76,61 @@ for (ss in 1:length(sim.seed)){ # tic() print(paste0("Parameter set: ", ss)) vparameters.temp <- sim.data.ls[[ss]] - vparameters.temp$NxDataPharm$pe <- round(vparameters.temp$NxDataPharm$pe / 2.69, 0) # ATTN: This is temporary, used to manually adjust pharmacy naloxone - # vparameters.temp$gw.m.2inact <- 0.0059 # ATTN: This is temporary, used to manually add monthly growth for transition into inactive (for MAT increase) - vparameters.temp$gw.m.2inact <- 0 # ATTN: This is temporary, used to manually add monthly growth for transition into inactive (for MAT increase) - vparameters.temp$nlx.adj <- 1.2 # ATTN: This is temporary, used to manually adjust naloxone adjustment term to 1 - vparameters.temp$mor_bl <- 0.0588 * 0.58 - vparameters.temp$mor_nx <- 0.00899 * 0.58 - vparameters.temp$OD_cess <- 0 - vparameters.temp$gw.fx <- 0.05 - vparameters.temp$multi.relap <- 3 - vparameters.temp$multi.fx <- 6 - # vparameters.temp$p.preb2inact.ini <- vparameters.temp$p.preb2inact - vparameters.temp$p.preb2inact.ini <- 0.005 - vparameters.temp[(names(vparameters.temp)=='p.preb2inact')] <- NULL # ATTN: This is temporary, used to fix a recent change to - vparameters.temp$p.il.lr2inact.ini <- vparameters.temp$p.il.lr2inact - vparameters.temp[(names(vparameters.temp)=='p.preb2inact')] <- NULL # ATTN: This is temporary, used to manually adjust naloxone adjustment term to 1 - vparameters.temp$p.il.hr2inact.ini <- vparameters.temp$p.il.hr2inact - vparameters.temp[(names(vparameters.temp)=='p.preb2inact')] <- NULL # ATTN: This is temporary, used to manually adjust naloxone adjustment term to 1 - vparameters.temp$Low2Priv <- 0.3 - vparameters.temp$High2Pub <- 0.1 - # vparameters.temp$p.inact2relap <- 0.004 - # vparameters.temp$overdose_probs = vparameters.temp$overdose_probs * 0.9 - - # status quo scenario 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 - sim_sq <- MicroSim(init_ppl, params = params, timesteps, agent_states, discount.rate, PT.out = FALSE, strategy = "SQ", seed = 2002) # run for status quo - + mx.od.death.1st[ss, "Status Quo"] <- sum(sim_sq$m.oddeath[(timesteps-35):(timesteps-24), ]) + mx.od.death.2nd[ss, "Status Quo"] <- sum(sim_sq$m.oddeath[(timesteps-23):(timesteps-12), ]) mx.od.death.last[ss, "Status Quo"] <- sum(sim_sq$m.oddeath[(timesteps-11):timesteps, ]) mx.od.death.totl[ss, "Status Quo"] <- sum(sim_sq$m.oddeath[(timesteps-35):timesteps, ]) + mx.od.death.wtns.1st[ss, "Status Quo"] <- sum(sim_sq$v.oddeath.w[(timesteps-35):(timesteps-24)]) + mx.od.death.wtns.2nd[ss, "Status Quo"] <- sum(sim_sq$v.oddeath.w[(timesteps-23):(timesteps-12)]) mx.od.death.wtns.last[ss, "Status Quo"] <- sum(sim_sq$v.oddeath.w[(timesteps-11):timesteps]) mx.od.death.wtns.totl[ss, "Status Quo"] <- sum(sim_sq$v.oddeath.w[(timesteps-35):timesteps]) mx.costs.totl[ss, "Status Quo"] <- sim_sq$total.cost array.od.death.rgn[ , "Status Quo", ss] <- colSums(sim_sq$m.oddeath[(timesteps-11):timesteps, ]) - array.oend.nlx.rgn[ , "Status Quo", ss] <- sim_sq$n.nlx.OEND.str + array.oend.nlx.rgn[ , "Status Quo", ss] <- as.vector(t(sim_sq$n.nlx.OEND.str)) for (jj in 2:length(scenario.name)){ vparameters.temp$expand.kits <- 10000 sim_pg <- MicroSim(init_ppl, params = vparameters.temp, timesteps, agent_states, discount.rate, PT.out = FALSE, strategy = scenario.name[jj], seed = sim.seed[ss]) # run for program scenario + mx.od.death.1st[ss, jj] <- sum(sim_pg$m.oddeath[(timesteps-35):(timesteps-24), ]) + mx.od.death.2nd[ss, jj] <- sum(sim_pg$m.oddeath[(timesteps-23):(timesteps-12), ]) mx.od.death.last[ss, jj] <- sum(sim_pg$m.oddeath[(timesteps-11):timesteps, ]) mx.od.death.totl[ss, jj] <- sum(sim_pg$m.oddeath[(timesteps-35):timesteps, ]) + mx.od.death.wtns.1st[ss, jj] <- sum(sim_pg$v.oddeath.w[(timesteps-35):(timesteps-24)]) + mx.od.death.wtns.2nd[ss, jj] <- sum(sim_pg$v.oddeath.w[(timesteps-23):(timesteps-12)]) mx.od.death.wtns.last[ss, jj] <- sum(sim_pg$v.oddeath.w[(timesteps-11):timesteps]) mx.od.death.wtns.totl[ss, jj] <- sum(sim_pg$v.oddeath.w[(timesteps-35):timesteps]) mx.costs.totl[ss, jj] <- sim_pg$total.cost array.od.death.rgn[ , jj, ss] <- colSums(sim_pg$m.oddeath[(timesteps-11):timesteps, ]) - array.oend.nlx.rgn[ , jj, ss] <- sim_pg$n.nlx.OEND.str + array.oend.nlx.rgn[ , jj, ss] <- as.vector(t(sim_pg$n.nlx.OEND.str)) } # toc() } +rgn.results <- data.frame(location = rep(v.region, 4), scenario = rep(scenario.name, each=num_regions), + N_Nx_mean = numeric(num_regions * 4), N_Nx_lower = numeric(num_regions * 4), N_Nx_upper = numeric(num_regions * 4), + N_ODDeath_mean = numeric(num_regions * 4), N_ODDeath_lower = numeric(num_regions * 4), N_ODDeath_upper = numeric(num_regions * 4), + population = rep(colSums(Demographic[,-c(1:3)]), 4)) -results.od.death.rgn <- apply(array.od.death.rgn, c(1,2), function(x) quantile(x, probs = c(0.025, 0.5, 0.975))) -results.oend.nlx.rgn <- apply(array.oend.nlx.rgn, c(1,2), function(x) quantile(x, probs = c(0.025, 0.5, 0.975))) -new.column.name <- as.vector(t(outer(dimnames(results.od.death.rgn)[[3]], dimnames(results.od.death.rgn)[[1]], paste, sep="_"))) -for(ii in 1:dim(results.od.death.rgn)[3]){ - if (ii == 1){ - mx.results.od.death.rgn <- t(results.od.death.rgn[,,ii]) - mx.results.oend.nlx.rgn <- t(results.oend.nlx.rgn[,,ii]) - } else { - mx.results.od.death.rgn <- cbind(mx.results.od.death.rgn, t(results.od.death.rgn[,,ii])) - mx.results.oend.nlx.rgn <- cbind(mx.results.oend.nlx.rgn, t(results.oend.nlx.rgn[,,ii])) - } +for (ii in 1:length(scenario.name)){ + rgn.results$N_Nx_mean[rgn.results$scenario== scenario.name[ii]] <- apply(array.oend.nlx.rgn, c(1,2), mean)[,scenario.name[ii]] + rgn.results$N_Nx_lower[rgn.results$scenario== scenario.name[ii]] <- apply(array.oend.nlx.rgn, c(1,2), function(x) quantile(x, probs = c(0.025)))[,scenario.name[ii]] + rgn.results$N_Nx_upper[rgn.results$scenario== scenario.name[ii]] <- apply(array.oend.nlx.rgn, c(1,2), function(x) quantile(x, probs = c(0.975)))[,scenario.name[ii]] + rgn.results$N_ODDeath_mean[rgn.results$scenario== scenario.name[ii]] <- apply(array.od.death.rgn, c(1,2), mean)[,scenario.name[ii]] + rgn.results$N_ODDeath_lower[rgn.results$scenario== scenario.name[ii]] <- apply(array.od.death.rgn, c(1,2), function(x) quantile(x, probs = c(0.025)))[,scenario.name[ii]] + rgn.results$N_ODDeath_upper[rgn.results$scenario== scenario.name[ii]] <- apply(array.od.death.rgn, c(1,2), function(x) quantile(x, probs = c(0.975)))[,scenario.name[ii]] } -colnames(mx.results.od.death.rgn) <- colnames(mx.results.oend.nlx.rgn) <- new.column.name -ppl_region <- colSums(Demographic[, -c(1:3)]) -mx.results.od.death.rgn <- cbind(mx.results.od.death.rgn, ppl_region) -mx.results.oend.nlx.rgn <- cbind(mx.results.oend.nlx.rgn, ppl_region) +detach("package:openxlsx", unload = TRUE) +library(xlsx) -write.xlsx(mx.od.death.last, file = ("Outputs/RI10K/ODdeaths_LastYear.xlsx"), row.names = F) -write.xlsx(mx.od.death.totl, file = ("Outputs/RI10K/ODdeaths_Total3Y.xlsx"), row.names = F) -write.xlsx(mx.od.death.wtns.last, file = ("Outputs/RI10K/WitnessedDeaths_LastYear.xlsx"), row.names = F) -write.xlsx(mx.od.death.wtns.totl, file = ("Outputs/RI10K/WitnessedDeaths_Total3Y.xlsx"), row.names = F) +write.xlsx(mx.od.death.1st, file = ("Outputs/RI10K/ODdeaths.xlsx"), sheetName = "1st year", row.names = F) +write.xlsx(mx.od.death.2nd, file = ("Outputs/RI10K/ODdeaths.xlsx"), sheetName = "2nd year", append = T, row.names = F) +write.xlsx(mx.od.death.last, file = ("Outputs/RI10K/ODdeaths.xlsx"), sheetName= "last year", append = T, row.names = F) +write.xlsx(mx.od.death.totl, file = ("Outputs/RI10K/ODdeaths.xlsx"), sheetName= "Total", append = T, row.names = F) +write.xlsx(mx.od.death.wtns.1st, file = ("Outputs/RI10K/WitnessedDeaths.xlsx"), sheetName = "1st year", row.names = F) +write.xlsx(mx.od.death.wtns.2nd, file = ("Outputs/RI10K/WitnessedDeaths.xlsx"), sheetName = "2nd year", append = T, row.names = F) +write.xlsx(mx.od.death.wtns.last, file = ("Outputs/RI10K/WitnessedDeaths.xlsx"), sheetName= "last year", append = T, row.names = F) +write.xlsx(mx.od.death.wtns.totl, file = ("Outputs/RI10K/WitnessedDeaths.xlsx"), sheetName= "Total", append = T, row.names = F) write.xlsx(mx.costs.totl, file = ("Outputs/RI10K/Total_costs.xlsx"), row.names = F) -write.xlsx(mx.results.od.death.rgn, file = ("Outputs/RI10K/ODdeaths_bytown.xlsx"), row.names = T) -write.xlsx(mx.results.oend.nlx.rgn, file = ("Outputs/RI10K/OENDnaloxone_bytown.xlsx"), row.names = T) +write.csv(rgn.results, file = ("Outputs/RI10K/Results_bytown.csv"), row.names = F) diff --git a/naloxone_availability.R b/naloxone_availability.R index faa6fa2..e9bb624 100644 --- a/naloxone_availability.R +++ b/naloxone_availability.R @@ -28,8 +28,7 @@ nlx.avail.algm <- function(n.nlx, crosstab_drug_resid, OD_loc_pub, nlx.adj, cap, } else if(t <= 48){ n.nlx$OEND <- n.nlx$OEND[-length(n.nlx$OEND)] n.nlx.total <- n.nlx$OEND + n.nlx$Pharm / eff.pharNlx - n.nlx.drug.resid <- as.vector(t(n.nlx.total)) * crosstab_drug_resid / rowSums(crosstab_drug_resid) - p.nlx.avail <- cap * (1 - exp(-(nlx.adj / cap) * n.nlx.drug.resid / rowSums(crosstab_drug_resid))) + p.nlx.avail <- cap * (1 - exp(-(nlx.adj / cap) * n.nlx.total / rowSums(crosstab_drug_resid))) return(p.nlx.avail) } else { tenk <- sum(n.nlx$OEND[length(n.nlx$OEND)]) @@ -39,16 +38,7 @@ nlx.avail.algm <- function(n.nlx, crosstab_drug_resid, OD_loc_pub, nlx.adj, cap, SSP_drug_resid <- crosstab_drug_resid[, c("il.hr", "NODU")]*rep(c(1, 0.133), each = nrow(crosstab_drug_resid)) n.nlx.drug.resid <- n.nlx.drug.resid.base n.nlx.drug.resid[, c("il.hr", "NODU")] <- tenk * SSP_drug_resid / sum(SSP_drug_resid) + n.nlx.drug.resid.base[, c("il.hr", "NODU")] - # } else if (strategy == "StreetOutreach_10K"){ - # StreetOutreach_pub_drug_resid <- pub_drug_resid[, c("il.lr", "il.hr", "NODU")] - # n.nlx.drug.resid.priv <- n.nlx.drug.resid.base.priv; n.nlx.drug.resid.pub <- n.nlx.drug.resid.base.pub - # n.nlx.drug.resid.pub[, c("il.lr", "il.hr", "NODU")] <- tenk * StreetOutreach_pub_drug_resid / sum(StreetOutreach_pub_drug_resid) + n.nlx.drug.resid.base.pub[, c("il.lr", "il.hr", "NODU")] } else if (strategy == "MailEvent_10K"){ - # MailEvent_priv_drug_resid <- priv_drug_resid - # MailEvent_pub_drug_resid <- pub_drug_resid - # n.nlx.drug.resid.priv <- n.nlx.drug.resid.base.priv; n.nlx.drug.resid.pub <- n.nlx.drug.resid.base.pub - # n.nlx.drug.resid.priv <- tenk * Low2Priv * MailEvent_priv_drug_resid / sum(MailEvent_priv_drug_resid) + n.nlx.drug.resid.base.priv - # n.nlx.drug.resid.pub <- tenk * (1-Low2Priv) * MailEvent_pub_drug_resid / sum(MailEvent_pub_drug_resid) + n.nlx.drug.resid.base.pub n.nlx.drug.resid <- n.nlx.drug.resid.base n.nlx.drug.resid <- tenk * crosstab_drug_resid / sum(crosstab_drug_resid) + n.nlx.drug.resid.base } else if (strategy == "Healthcare_10K"){