-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathpart_3.R
118 lines (103 loc) · 4.22 KB
/
part_3.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
##########################################################################
# Run some more nuanced simulations using a wrapper function for nimue::run, which allows us to set some more flexible options:
# - income setting (or change demography/contact patterns/population size)
# - transmission over time
# - vaccine coverage
# - vaccine dose strategy (age groups targeted and prioritisation)
# - vaccines delivered per day
# - vaccine efficacy
# - whether vaccine targets infection, disease, or both
# - time period over which simulation is run
# - hospital capacity
##########################################################################
# Now we will look at an example with some different vaccine options
# Load the functions script
source("R/functions.R")
# set some parameters
R0 <- 2.5
Rt1 <- 0.9
Rt2 <- 3
reduction1 <- 1-Rt1/R0
reduction2 <- 1-Rt2/R0
timing1 <- 50
timing2 <- 350
income_group <- "HIC"
duration_R <- 365
coverage <- c(0.8,0.6, 0)
vaccine_coverage_mat <- "Elderly"
vaccine_period <- 30
vaccine_start <- 365
# set up a dataframe of scenarios
scenarios <- expand_grid(R0 = R0,
reduction1 = reduction1,
reduction2 = reduction2,
coverage = coverage,
income_group = income_group,
duration_R = duration_R,
vaccine_period = vaccine_period,
vaccine_start = vaccine_start,
vaccine_coverage_mat = vaccine_coverage_mat,
timing1 = timing1,
timing2 = timing2)
# run the model
system.time({out <- future_pmap(scenarios, run_scenario, .progress = TRUE)})
# format the output
df <- bind_cols(scenarios, bind_rows(out)) %>%
unnest(output) %>%
filter(compartment %in% c("deaths", "infections"))
# create plot
ggplot(data = df, aes(x = t, y = value, col = factor(coverage))) +
geom_line(size = 0.8) +
facet_wrap(~compartment, scales = "free") +
theme_bw() +
scale_color_viridis_d(end = 0.8)
# what about total events (deaths or infections) over the time period?
df_summary <- df %>%
group_by(coverage, compartment) %>%
summarise(value = sum(value, na.rm = T))
ggplot(data = df_summary, aes(x = factor(coverage), y = value, fill = factor(coverage))) +
geom_bar(stat = "identity") +
facet_wrap(~compartment, scales = "free") +
theme_bw() +
scale_fill_viridis_d(end = 0.8) +
labs(x = "coverage", y = "total events")
###########################################################################################
# what about impact of changing vaccine efficacy?
coverage <- 0.6
efficacy <- seq(0.3, 1, by = 0.05)
scenarios <- expand_grid(R0 = R0,
reduction1 = reduction1,
reduction2 = reduction2,
coverage = coverage,
efficacy = efficacy,
income_group = income_group,
duration_R = duration_R,
vaccine_period = vaccine_period,
vaccine_start = vaccine_start,
vaccine_coverage_mat = vaccine_coverage_mat,
timing1 = timing1,
timing2 = timing2)
# run the model
system.time({out <- future_pmap(scenarios, run_scenario, .progress = TRUE)})
# format the output
df <- bind_cols(scenarios, bind_rows(out)) %>%
unnest(output) %>%
filter(compartment %in% c("deaths", "infections"))
# create plot
ggplot(data = df, aes(x = t, y = value, col = factor(efficacy))) +
geom_line(size = 0.8) +
facet_wrap(~compartment, scales = "free") +
theme_bw() +
scale_color_viridis_d(end = 0.8)
# what about total events (deaths or infections) over the time period?
df_summary <- df %>%
group_by(efficacy, compartment) %>%
summarise(value = sum(value, na.rm = T))
ggplot(data = df_summary, aes(x = factor(efficacy), y = value, fill = factor(efficacy))) +
geom_bar(stat = "identity") +
facet_wrap(~compartment, scales = "free") +
theme_bw() +
scale_fill_viridis_d(end = 0.8) +
labs(x = "coverage", y = "total events")
###########################################################################################
# what about impact in different income settings?