-
Notifications
You must be signed in to change notification settings - Fork 1
/
README.Rmd
277 lines (205 loc) · 9.12 KB
/
README.Rmd
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
*simex*: a disease modelling tool for decision makers
---------------------------------------------------------------
[](https://github.com/epiforecasts/EpiNow2/blob/main/LICENSE.md/) [](https://lifecycle.r-lib.org/articles/stages.html#maturing)
*simex* is an R package for simulating the spread of infectious diseases, as well as interventions such as social distancing measures, isolation and vaccination. It uses an age-structured SEIR compartmental model with country-specific age demographics and contact rates. The simplest way to interact with the tool is to use the [Shiny App](https://portal.who.int/eios-colab/rconnect/simex).
```{r, echo = FALSE}
knitr::opts_chunk$set(fig.cap = "",
fig.align = "center",
dpi = 200,
out.width = "75%")
```
Installation
-------------
To install the development version from github:
```{r, eval = FALSE}
remotes::install_github("WHO-Collaboratory/simex")
```
Load the package using:
```{r}
library("simex")
```
Running *simex*
-------------
### Shiny App
You can use *simex* interactively by launching the shiny app locally with
`simex::run_shiny()`. If you want to use it programmatically, follow the steps
below.
### Parameters and settings
Most settings are specified via the `get_parameters` function. The arguments and
their default values are given below:
```{r, echo = FALSE}
simex:::get_arg_table("get_parameters", type = "kable")
```
### Running default settings
To run the model using default settings, specify a parameter object `par` and
feed this into the `run_model` function.
```{r, echo = TRUE}
## set parameters using defaults
pars <- get_parameters()
## run model using default parameters
output <- run_model(pars)
## look at output
print(output)
```
### Visualising outputs
To visualise the results, use the generic `plot` function defined for the
`simex` class. Below, we first visualise prevalence and then incidence,
specified using the `what` argument.
```{r, echo = TRUE}
## visualise prevalence
plot(output, what = "prevalence")
## visualise incidence
plot(output, what = "incidence")
```
Hospital capacity can be displayed by toggling the `show_hosp_capacity`
argument.
```{r, echo = TRUE}
## visualise prevalence with hospital capaciy
plot(output, what = "prevalence", show_hosp_capacity = TRUE)
```
### Summarising outputs
To get a summary of the model state at a given point in time, use the `summary`
function with the optional `day` argument. In the example below, we display the
summary statistics for day 150:
```{r, eval = FALSE}
## get summary at day 150
summary(output, day = 150)
```
```{r, echo = FALSE}
## kable for Rmd
knitr::kable(summary(output, day = 150))
```
### Accessing outputs
The `output` object is of class `simex` and is a list containing four items:
* `prevalence` is an array with 3 dimensions: time (365 days) x age (16
categories) x state (12 compartments). The compartments are S (susceptible), E
(exposed, pre-symptomatic), C (symptomatic in community), H (sympomatic in
hospital), R (recovered), D (dead). Each comparment is split into unvaccinated
(with subscript _u) and vaccinated (with subscript _v). Each value in the
array represents the proportion of the population in that given age-disease
compartment on a given day.
* `deltas` is an array with the same dimensions as `prevalence`. Each value
represents the change in a given age-disease compartment from the previous
day.
* `incidence` is an array with the same dimensions as `prevalence`. Each value
represents the new additions to a given age-disease comparment on a given
day. For the exposed category E, this represents incidence of new infections,
for the hospitalised category H this represents new hospital admissions, and
for the dead category D this represents new deaths.
* `pars` contains the parameter set initially fed into `run_model`.
Accessing the data is easiest using array indexing. Remember, the dimensions are
time, age and compartment. For example, accessing the prevalence on the 250th
day is done with `output$prevalence[250,,]`:
```{r, echo = FALSE}
## extract prevalence in percent on day 250
round(100*output$prevalence[250,,], 1)
```
Accessing the prevalence of the 1st age compartment (0-4) and 1st infectious
compartment (unvaccinated susceptible) for days 130 to 135 is done with
`output$prevalence[130:135,1,1]`:
```{r, echo = FALSE}
## extract in percent on day 130-135
round(100*output$prevalence[130:135,1,1], 1)
```
The outputs can also be extracted in `tibble` form using the `extract` function,
once again using the `what` argument to specify whether prevalence or incidence
is extracted.
```{r, echo = TRUE}
## extract prevalence
extract(output, what = "prevalence")
```
If we want to filter this list for a sequence of days, we can then do basic
dataframe manipulation:
```{r, echo = TRUE}
## define start and end days
days_from <- 10
days_to <- 20
## extract prevalence
df <- extract(output, what = "prevalence")
df <- df[df$day %in% seq(days_from, days_to),]
df
```
### Modelling a single intervention
We use vaccination as an example intervention. Referencing the table above, we
can see that vaccination rate is specified using the `vax_rate` argument and
set it to 0.5% of the population per day.
```{r, echo = TRUE}
## define vaccination rate
pars <- get_parameters(vax_rate = 0.005)
## run model
output <- run_model(pars)
## visualise prevalence
plot(output)
```
We can see that the daily increase in number of vaccinated individuals, as
well as the impact on infection and disease severity.
### Specifying an initial state
The default initial state begins with an entirely susceptible, unvaccinated
population and a single infection. We can also specify a custom initial state by
passing a matrix specifying the number of individuals in each age-infection
compartment. In the example below, we extract the initial state from the
previous run, and modify the compartments so half the susceptible population is
assigned to the vaccinated compartment.
```{r, echo = TRUE}
## extract starting point from previous run
state <- output$prevalence[1,,]
## assign half of susceptibles to vaccinated
state[,"S_u"] <- state[,"S_v"] <- state[,"S_u"]/2
## run model
output <- run_model(pars, init_state = state)
## visualise prevalence
plot(output)
```
We can see that the simulation now begins with a 50% vaccinated population,
significantly reducing the size and impact of the pandemic.
### Timed interventions, multiple interventions
Sometime we only want to introduce interventions at a certain time, or we want
to model multiple, successive interventions. We can do this easily by generating
a named list of parameters, where the name gives the time that parameter set
should be used from. In the example below, we introduce isolation measures with
an adherence of 50% on the 125th day.
```{r, echo = TRUE}
## introduce isolation on the 125th day
parlist <- list(
"1" = get_parameters(),
"125" = get_parameters(isolation_adherence = 0.5)
)
## run model
output <- run_model(parlist)
## visualise prevalence
plot(output)
```
Comparing this figure with the first model run with no interventions, we can see
the proportion of deaths drops from about 0.7% to 0.4%; a reduction in deaths of
more than 40%!
### Comparing scenarios
It is useful to directly compare different scenarios visually. The
`vis_comarison` function does exactly this, and accepts a named list of `simex`
objects. In the example below, we generate a named list of lists that compares
the default scenario (no interventions) with the a scenario where isolation is
introduced on the 125th day.
```{r, echo = TRUE}
## define two scenarios, one without intervention and one with isolation
parlists <- list(
"No intervention" = get_parameters(),
"Isolation on day 125" = list(
"1" = get_parameters(),
"125" = get_parameters(isolation_adherence = 0.5)
)
)
## run model across set of parameter lists
outputs <- lapply(parlists, run_model)
## compare scenarios
vis_comparison(outputs)
```
### Contributors
------------
- [Finlay Campbell](https://github.com/finlaycampbell) ([email protected])
- Prabasaj Paul ([email protected])
**Maintainer:** Finlay Campbell
### Licensing
The simex software is made available by Collaboratory (2024) under an [MIT license](LICENSE.md) ([citation file](CITATION.cff)).
The demographic data used in this software is made available by United Nations (2024) under a [Creative Commons license CC BY 3.0 IGO](http://creativecommons.org/licenses/by/3.0/igo/). Source:
> United Nations, Department of Economic and Social Affairs, Population Division (2024). [World Population Prospects 2024, Online Edition](https://population.un.org/wpp/).
The contact data used in this software is made available by Prem et al. (2017) under a [Creative Commons license CC BY 4.0](https://creativecommons.org/licenses/by/4.0/). Source:
> Kiesha Prem, Alex R. Cook, Mark Jit, *Projecting social contact matrices in 152 countries using contact surveys and demographic data*, PLoS Comp. Biol. (2017), https://doi.org/10.1371/journal.pcbi.1005697.