-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathrun_dynamice_201910gavi.R
284 lines (242 loc) · 11.6 KB
/
run_dynamice_201910gavi.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
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
# run_dynamice.R
# run it from the source directory
# DynaMICE (Dynamic Measles Immunisation Calculation Engine)
# Measles vaccine impact model
# Estimate health impact (cases, deaths, dalys) of measles vaccination for a
# given set of countries.
# remove all objects from workspace
rm (list = ls ())
# ------------------------------------------------------------------------------
# load libraries (delete when package is ready)
library (tictoc)
library (data.table)
library (stringr)
library (ggplot2)
library (scales)
library (foreach)
library (iterators)
library (parallel)
library (doParallel)
library (countrycode)
library (ggpubr)
library (lhs)
library (truncnorm)
# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
# load data (delete when package is ready)
load(file = "data/data_pop.rda")
load(file = "data/data_cfr_portnoy.rda")
load(file = "data/data_cfr_wolfson.rda")
load(file = "data/data_contact_polymod.rda")
load(file = "data/data_contact_syn.rda")
load(file = "data/data_r0.rda")
load(file = "data/data_timeliness.rda")
load(file = "data/data_lexp_remain.rda")
load(file = "data/data_template.rda")
# load(file = "data/data_cfr.rda")
#load(file = "data/data_psa200.rda")
# load(file = "data/data_contact_uk.rda")
# load(file = "data/data_lexp0.rda")
# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
# include functions (delete when package is ready)
source ("R/logs.R")
source ("R/utils.R")
source ("R/functions.R")
source ("R/dxplot.R")
# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
# main program -- start
# ------------------------------------------------------------------------------
# start time
print (Sys.time ())
# set seed for random number generator
set.seed (1)
# ------------------------------------------------------------------------------
# set values for global variables
var <- list (
# vaccine coverage
vaccine_coverage_folder = "vaccine_coverage/",
coverage_prefix = "coverage",
touchstone = "_201910gavi-5_",
antigen = "measles-",
vaccine_coverage_subfolder = "scenarios/",
# disease burden
central_burden_estimate_folder = "central_burden_estimate/",
stochastic_burden_estimate_folder = "stochastic_burden_estimate/",
# diagnostic plots folder
plot_folder = "plots/",
# modelling group name
group_name = "LSHTM-Jit-",
# countries - specify iso3 codes to analyse only these countries
# or set it to "all" to analyse all included countries
countries = c("BGD","PAK"), # debug -- c( "ETH") / "all"
cluster_cores = 3, # number of cores
psa = 0 # psa runs; 0 for single central run
)
# for central run: set number of runs to 0 (var$psa = 0)
# for stochastic runs: set number of runs to greater than 0 (eg: var$psa = 200)
# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
# scenarios
scenarios <- c("campaign-only-bestcase", # 1 SIAs only
"mcv2-bestcase", # 2 MCV1 & MCV2
"campaign-bestcase", # 3 MCV1 & MCV2 and SIAs
"mcv1-bestcase", # 4 MCV1 only
"campaign-only-default", # 5 SIAs only
"mcv2-default", # 6 MCV1 & MCV2
"campaign-default", # 7 MCV1 & MCV2 and SIAs
"mcv1-default", # 8 MCV1 only
"no-vaccination", # 9 no vaccination (set vaccination and using_sia to 0)
"stop" # 10 MCV1 & MCV2 and SIAs
)
# ------------------------------------------------------------------------------
# specify scenarios to run
first_scenario <- 1
last_scenario <- length (scenarios)
# ------------------------------------------------------------------------------
# set SIAs and vaccination parameters for each scenario to minimize errors for running
set_sia <- c (1, 0, 1, 0, 1, 0, 1, 0, 0, 1)
set_vaccination <- c (0, 2, 2, 1, 0, 2, 2, 1, 0, 2)
# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
# Create vaccine coverage file (0% coverage) for no vaccination scenario using
# another vaccination scenario. This is done because VIMC vaccine coverage file
# for no vaccination scenario is empty
create_no_vaccination_coverage_file (
no_vaccination_coverage_file = "vaccine_coverage/coverage_201910gavi-5_measles-no-vaccination.csv",
vaccination_coverage_file = "vaccine_coverage/coverage_201910gavi-5_measles-mcv1-default.csv"
)
# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
# Create vaccine coverage file for campaign only vaccination scenario using
# (routine + campaign) vaccination scenario.
# Set routine coverage to zero. This is done because routine coverage values are
# needed even if they are only zero to run campaign only vaccination scenario.
# scenario: campaign-only-bestcase
create_campaign_vaccination_coverage_file (
campaign_only_vaccination_coverage_file = "vaccine_coverage/coverage_201910gavi-5_measles-campaign-only-bestcase.csv",
routine_campaign_vaccination_coverage_file = "vaccine_coverage/coverage_201910gavi-5_measles-campaign-bestcase.csv"
)
# scenario: campaign-only-default
create_campaign_vaccination_coverage_file (
campaign_only_vaccination_coverage_file = "vaccine_coverage/coverage_201910gavi-5_measles-campaign-only-default.csv",
routine_campaign_vaccination_coverage_file = "vaccine_coverage/coverage_201910gavi-5_measles-campaign-default.csv"
)
# ------------------------------------------------------------------------------
# Generate 2 vaccine coverage files per scenario for routine and SIA from
# VIMC vaccine coverage files
for (index in 1:length(scenarios)) {
create_vaccine_coverage_routine_sia (
vaccine_coverage_folder = var$vaccine_coverage_folder,
vaccine_coverage_subfolder = var$vaccine_coverage_subfolder,
coverage_prefix = var$coverage_prefix,
touchstone = var$touchstone,
antigen = var$antigen,
scenario_name = scenarios [index],
rev_cov = TRUE
)
}
# ------------------------------------------------------------------------------
# # create remaining life expectancy file for each year across all age intervals
# create_life_expectancy_remaining_full ()
# create latin hyper cube sample of input parameters for
# probabilistic sensitivity analysis
if (var$psa > 0) {
psadat <- create_PSA_data (psa = var$psa,
seed_state = 1,
psadat_filename = "psa_variables.csv")
}
# burden_estimate_folder (central or stochastic)
if (var$psa == 0) {
burden_estimate_folder <- var$central_burden_estimate_folder
} else {
burden_estimate_folder <- var$stochastic_burden_estimate_folder
}
# loop through scenarios first_scenario:last_scenario, c(1,3,5,7,9,10)
for (index in 6:7) {
# set scenario name
scenario_name <- scenarios [index]
print (scenario_name)
tictoc::tic ()
# set scenario number in order to save it in the right folder
# scenarioXX, XX = 01, 02, ... ,10, 11, 12, ...
# this is set to 10 characters in the fortran model code
scenario_number <- sprintf ("scenario%02d", index)
# ----------------------------------------------------------------------------
# estimate cases
#
# run scenario -- get burden estimates -- primarily cases
# return burden estimate file name where estimates are saved
burden_estimate_file <- runScenario (
vaccine_coverage_folder = var$vaccine_coverage_folder,
vaccine_coverage_subfolder = var$vaccine_coverage_subfolder,
coverage_prefix = var$coverage_prefix,
antigen = var$antigen,
scenario_name = scenario_name,
save_scenario = scenario_number,
burden_estimate_folder = burden_estimate_folder,
group_name = var$group_name,
countries = var$countries,
cluster_cores = var$cluster_cores,
psa = var$psa,
vaccination = set_vaccination [index],
using_sia = set_sia [index],
measles_model = "vaccine2019_sia_adjBeta_mcv2.exe",
debug_model = FALSE,
contact_mat = "polymod",
step_ve = FALSE
)
# ----------------------------------------------------------------------------
# ----------------------------------------------------------------------------
# estimate deaths and DALYs
#
# apply CFR (case fatality rates) to estimate deaths -- Wolfson
# save results in corresponding cfr_option subfolder
# append cfr_option to results file
# uncomment below to run with cfr_option = "Wolfson"
estimate_deaths_dalys (cfr_option = "Wolfson",
burden_estimate_file = burden_estimate_file,
burden_estimate_folder = burden_estimate_folder,
psa = var$psa)
# apply CFR (case fatality rates) to estimate deaths -- Portnoy
# save results in corresponding cfr_option subfolder
# append cfr_option to results file
estimate_deaths_dalys (cfr_option = "Portnoy",
burden_estimate_file = burden_estimate_file,
burden_estimate_folder = burden_estimate_folder,
vimc_scenario = scenario_name,
portnoy_scenario = "s6", # portnoy scenario 6
psa = var$psa)
# ----------------------------------------------------------------------------
tictoc::toc ()
} # end of loop -- for (scenario in scenarios)
# base scenario for comparison
base_scenario <- "no-vaccination"
# ------------------------------------------------------------------------------
# diagnostic plots of vaccine coverage and burden estimates (cases, deaths, dalys)
if (var$psa == 0) {
diagnostic_plots (
vaccine_coverage_folder = var$vaccine_coverage_folder,
coverage_prefix = var$coverage_prefix,
touchstone = var$touchstone,
antigen = var$antigen,
scenarios = scenarios [first_scenario:last_scenario],
base_scenario = base_scenario,
burden_estimate_folder = var$central_burden_estimate_folder,
plot_folder = var$plot_folder,
group_name = var$group_name,
countries = var$countries,
cfr_options = c("Wolfson", "Portnoy"),
psa = var$psa,
start_year = 2000,
end_year = 2050,
compare_plots = FALSE
)
}
# ------------------------------------------------------------------------------
# end time
print (Sys.time ())
# ------------------------------------------------------------------------------
# main program -- end
# ------------------------------------------------------------------------------