-
Notifications
You must be signed in to change notification settings - Fork 1
/
ms7.2.1.multivar_countryLTs_discount_RIG_ICER.R
166 lines (134 loc) · 7.95 KB
/
ms7.2.1.multivar_countryLTs_discount_RIG_ICER.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
################################################################################
# 7.2.1 Run Multivar decision tree model #
# OUTPUTS USED BY GAVI AND IN PAPER FOR BURDEN ESTIMATES #
# run ICER scenarios (discount) i.e. baseline improved PEP, RIG & dog vax
# save outputs to folder: countryLTs_discount
################################################################################
#' * Life Tables - Country specific *
#' * Discounting - 3% *
#' * PEP cost - $5 (default) *
#' * RIG cost - $45 (default) *
#' * Intro grant - $100k *
#' * Scenarios - a1, a3_1, a4, a2, a5_1, a5_2 *
#' * Run count - 1000 *
rm(list=ls())
# Load in packages
library(gdata)
library(rlang)
library(reshape2)
library(ggplot2)
library(tools)
library(triangle)
library(plyr)
library(dplyr)
library(Hmisc)
library(tidyr)
# Load in functions
source("R/YLL.R") # Calculate YLL given life tables and rabies age distribution
source("R/PEP.R") # Vial use under different regimens and throughput
source("R/prob_rabies.R") # Probability of developing rabies - sensitivity analysis
source("R/decision_tree_sensitivity_by_year.R") # Sensitivity analysis
source("R/decision_tree_multivariate_analysis_by_year_v2.R") # Multivariate sensitivity analysis
source("R/scenario_params.R") # Parameters and functions for gavi support and phasing
source("R/multivar_output_summary_Github.R") # MEAN
source("R/multivariate_plot_summarise_data_Github.R")
# Set folder name for output
folder_name <- "countryLTs_discount"
######################
# 1. Setup variables #
######################
rabies = read.csv("data/baseline_incidence_Gavi_final.csv")
data <- read.csv("output/gavi_output_data.csv") # Load gavi-prepared data
params <- read.csv("output/bio_data.csv") # parameters i.e. rabies transmission, prevention given incomplete PEP
vacc <- read.csv("data/vaccine_use.csv") # PEP scenarios - clinic throughput, regimen, completeness, vials, clinic visits:
dogs <- read.csv(file="output/dogs_pop_traj.csv", stringsAsFactors = FALSE) # dog pop 2018-2070 created in 6.elimination_traj.R
elimination_traj <- read.csv(file="output/rabies_traj.csv") # by year of global business plan
y1 = "2020"; yN = "2035"
pop = data[,grep(y1, names(data)):grep(yN, names(data))] # needs this format to combine with elimination trajectories!
# Set time horizon: from 2020 to 2035
hrz=length(2020:2035)
# Load in DALYs - disability weightings and lifetables
DALYrabies_input <- read.csv("data/DALY_params_rabies.csv") # from Knobel et al. 2005
# SPECIFIC PARAMETERS
# Life table
country_LE <- read.csv("data/lifetables_bycountry.csv")
country_LE <- country_LE[-which(country_LE$age_from == 100),]
LE2020 <- country_LE[which(country_LE$year == 2020),] # Use 2020 age distributions throughout!
# Set discounting rate
discount = 0.03
# Set prices (USD)
gavi_intro_grant <- 100000 # Intro grant
gavi_vaccine_price <- 5 # vaccine cost per vial
gavi_RIG_price <- 45 # ERIG cost per vial
################
# 2. Run model #
################
# Set number of runs
n = 1000 # 1.5 hrs per run, so ~8 hrs
# SQ - Paper S1
scenario_a1 <- multivariate_analysis(ndraw=n, horizon=hrz, GAVI_status="none", DogVax_TF=F, VaxRegimen="Updated TRC",
DALYrabies=DALYrabies_input, LE=LE2020, RIG_status="none", discount=discount, breaks="5yr", IBCM=FALSE)
# Improved PEP - Paper SC2 (base)
scenario_a3_1 <- multivariate_analysis(ndraw=n, horizon=hrz, GAVI_status="base", DogVax_TF=F, VaxRegimen="Updated TRC",
DALYrabies=DALYrabies_input, LE=LE2020, RIG_status="none", discount=discount, breaks="5yr", IBCM=FALSE)
# Improved PEP + RIG - Paper SC3
scenario_a4 <- multivariate_analysis(ndraw=n, horizon=hrz, GAVI_status="base", DogVax_TF=F, VaxRegimen="Updated TRC",
DALYrabies=DALYrabies_input, LE=LE2020, RIG_status="high risk", discount=discount, breaks="5yr", IBCM=FALSE)
# Dog vacc SQ - Paper SC4a
scenario_a2 <- multivariate_analysis(ndraw=n, horizon=hrz, GAVI_status="none", DogVax_TF=T, VaxRegimen="Updated TRC",
DALYrabies=DALYrabies_input, LE=LE2020, RIG_status="none", discount=discount, breaks="5yr", IBCM=FALSE)
# Dog vacc + PEP access - Paper SC4b
scenario_a5_1 <- multivariate_analysis(ndraw=n, horizon=hrz, GAVI_status="base", DogVax_TF=T, VaxRegimen="Updated TRC",
DALYrabies=DALYrabies_input, LE=LE2020, RIG_status="none", discount=discount, breaks="5yr", IBCM=FALSE)
# Dog vacc + PEP access + IBCM - Paper SC4c
scenario_a5_2 <- multivariate_analysis(ndraw=n, horizon=hrz, GAVI_status="base", DogVax_TF=T, VaxRegimen="Updated TRC",
DALYrabies=DALYrabies_input, LE=LE2020, RIG_status="none", discount=discount, breaks="5yr", IBCM=TRUE)
###########################################
# 3. Bind outputs into a single dataframe #
###########################################
# Append all results into a dataframe
out <- rbind.data.frame(
cbind.data.frame(scenario_a1, scenario="a1"),
cbind.data.frame(scenario_a3_1, scenario="a3_1"),
cbind.data.frame(scenario_a4, scenario="a4"),
cbind.data.frame(scenario_a2, scenario="a2"),
cbind.data.frame(scenario_a5_1, scenario="a5_1"),
cbind.data.frame(scenario_a5_2, scenario="a5_2"))
dim(out)
table(out$scenario)
countries <- unique(out$country)
scenarios <- unique(out$scenario)
yrs <- unique(out$year)
# INCLUDE GAVI ELIGIBILITY
gavi_info <- read.csv("output/gavi_output_data.csv", stringsAsFactors=FALSE)
out <- merge(out, data.frame(country=gavi_info$country, gavi_2018=gavi_info$gavi_2018), by="country", all.x=TRUE)
# CE outputs
out$cost_per_death_averted <- out$total_cost/out$total_deaths_averted
out$cost_per_YLL_averted <- out$total_cost/out$total_YLL_averted
out$deaths_averted_per_100k_vaccinated <- out$total_deaths_averted/out$vaccinated/100000
# Summarize by iteration over time horizon
out_horizon = country_horizon_iter(out)
######################################
# 4a. Create summary outputs #
######################################
# Country, cluster, & global by year
country_summary_yr = multivar_country_summary(out, year = TRUE)
cluster_summary_yr = multivar_summary(country_summary_yr, year=TRUE, setting ="cluster")
global_summary_yr = multivar_summary(country_summary_yr, year=TRUE, setting="global")
gavi2018_summary_yr = multivar_summary(country_summary_yr[which(country_summary_yr$gavi_2018==TRUE),], year=TRUE, setting="global")
write.csv(country_summary_yr, paste("output/", folder_name, "/country_stats_ICER.csv", sep=""), row.names=FALSE)
write.csv(cluster_summary_yr, paste("output/", folder_name, "/cluster_stats_ICER.csv", sep=""), row.names=FALSE)
write.csv(global_summary_yr, paste("output/", folder_name, "/global_stats_ICER.csv", sep=""), row.names=FALSE)
write.csv(gavi2018_summary_yr, paste("output/", folder_name, "/gavi2018_stats_ICER.csv", sep=""), row.names=FALSE)
################################################
# 4b. Create summary outputs over time horizon #
################################################
# Country, cluster, & global over time horizon
country_summary_horizon = multivar_country_summary(out_horizon, year = FALSE)
cluster_summary_horizon = multivar_summary(country_summary_horizon, year=FALSE, setting ="cluster")
global_summary_horizon = multivar_summary(country_summary_horizon, year=FALSE, setting="global")
gavi2018_summary_horizon = multivar_summary(country_summary_horizon[which(country_summary_horizon$gavi_2018==TRUE),], year=FALSE, setting="global")
write.csv(country_summary_horizon, paste("output/", folder_name, "/country_stats_horizon_ICER.csv", sep=""), row.names=FALSE)
write.csv(cluster_summary_horizon, paste("output/", folder_name, "/cluster_stats_horizon_ICER.csv", sep=""), row.names=FALSE)
write.csv(global_summary_horizon, paste("output/", folder_name, "/global_stats_horizon_ICER.csv", sep=""), row.names=FALSE)
write.csv(gavi2018_summary_horizon, paste("output/", folder_name, "/gavi2018_stats_horizon_ICER.csv", sep=""), row.names=FALSE)