-
Notifications
You must be signed in to change notification settings - Fork 1
/
part1.qmd
191 lines (137 loc) · 7.89 KB
/
part1.qmd
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
---
title: "Part 1: Basic Modeling"
---
**`epiworldR`** is an R package that provides a fast (C++ backend) and highly-customizable framework for building network-based transmission/diffusion agent-based models \[ABM\]. Some key features of **`epiworldR`** are the ability to construct multi-disease models (e.g., models of competing multi-pathogens/multi-rumor,) design mutating pathogens, architect population-level interventions, and build models with an arbitrary number of compartments/states (beyond SIR/SEIR.)[^1]
[^1]: This feature is currently under development. The [repository of epiworld](https://github.com/UofUEpiBio/epiworld){target="_blank"} contains a branch with this feature.
Let's start right away with an example!
## Motivating Example
```{r}
#| label: motivating-example
# Motivating example
library(epiworldR)
# Create a model
model <- ModelSEIRCONN(
name = "Monkeypox",
n = 50000,
prevalence = 0.0001,
contact_rate = 4,
incubation_days = 7,
transmission_rate = 0.5,
recovery_rate = 1/7
)
# Changing contact rate for Isolation and TV advertisement
isolation_day_10 <- globalaction_set_params("Contact rate", 2, day = 10)
advertisement_day_20 <- globalaction_set_params("Contact rate", 1.5, day = 20)
# Adding global actions to model
add_global_action(model, isolation_day_10)
add_global_action(model, advertisement_day_20)
# Running and printing model summary
run(model, ndays = 60, seed = 1912)
summary(model)
plot(model)
```
```{r}
#| label: run-model
#| eval: true
# Run the model
run(model, ndays = 100, seed = 1912)
model
```
```{r}
#| label: plot-model
# Plot the model
op <- par(mfrow = c(2,2))
plot_incidence(model)
abline(v = 20, col = "steelblue", lwd = 2, lty = 2)
plot_reproductive_number(model)
abline(v = 20, col = "steelblue", lwd = 2, lty = 2)
plot_generation_time(model)
abline(v = 20, col = "steelblue", lwd = 2, lty = 2)
par(op)
```
# General structure of epiworldR
epiworldR is a highly flexible framework for building agent-based-models. Although we do emphasize models of disease transmission, the framework is general enough to build models of any type of diffusion process. The general structure of epiworldR is as follows:
![A brief overview of epiworldR](figs/overview-of-epiworldr.drawio.svg){width="100%"}
## Example 1: Simulating a SIR model
#### Setup and running the model
This example implements the following scenario:
- The disease name is specified (COVID-19),
- 50,000 agents are initialized,
- the disease prevalence of 0.0001 is declared,
- each agent will contact two others (contact_rate),
- the transmission rate of the disease for any given agent is 0.5, and
- the recovery rate is set to $\frac{1}{3}$.
To create this model on **`epiworldR`**, we use the `ModelSIRCONN()` function. From here, the example will take you through the basic features of **`epiworldR`**.
```{r sirconn-setup}
library(epiworldR)
model_sirconn <- ModelSIRCONN(
name = "COVID-19",
n = 50000,
prevalence = 0.0001,
contact_rate = 2,
transmission_rate = 0.5,
recovery_rate = 1/3
)
```
Printing the model shows us some information. First, the name of the model, population size, number of entities (think of these as public spaces in which agents can make social contact with one another), the duration in days, number of variants, amount of time the last replicate took to run (last run elapsed t), and rewiring status (on or off). The next piece of information you will see is a list of the viruses used in the model. In this case, COVID-19 was the only disease used. Note that **`epiworldR`** has the capability to include more than one virus in a model. Tool(s) lists any tools that agents have to fight the virus. Examples ofthis may include masking, vaccines, social distancing, etc. In this model, no tools are specified. Lastly, the model parameters are listed, which originate from the parameters specified in the model.
```{r}
model_sirconn
```
To execute the model, use the run function with the SIR model object, number of simulation days, and an optional seed for reproducibility. Next, print out the results from the simulated model using model_sir.
```{r}
run(model_sirconn, ndays = 50, seed = 1912)
summary(model_sirconn)
```
There are two additional sections in the model summary after running the model object summary, the first being the distribution of the population at time 50. This section describes the flow of agents from each state (SIR) after 50 days. The counts for these states will of course, change based on model parameters or simulation run-time. The transmission probabilities section outputs a 3x3 matrix that describes the probability of moving from one state to another. Notice in all cases, there is a probability of 0 to skip states. In other words, it is impossible for an agent to move from the susceptible state to the recovered state; that agent must pass through the infected state in order to then progress to the recovered state. The same logic applies with moving backwards; an agent cannot become susceptible again after being infected.
#### Extracting Simulation Data
After running the **`epiworldR`** model, below is a list of all the functions that can be called using the **`epiworldR`** model object. To demonstrate, start with the basic plot and get_hist_total functions.
```{r showing-methods}
methods(class = "epiworld_model")
```
#### Visualization
```{r plotting}
plot(model_sirconn)
```
As evident from the above plot, the SIR model constructed from **`epiworldR`** displays the changes in susceptible, infected, and recovered case counts over time (days). Notice after a certain amount of time, the curves flatten. Below, a table representation of the above plot is printed, complete with each state within the SIR model, date, and agent counts.
::: {.callout-note collapse="true"}
## Connected vs Non-connected Models
The above example uses a SIR connected model (`ModelSIRCONN()`), meaning that all agents in the model are connected with each other. When using a non-connected model (ex.`ModelSEIR()`, `ModelSIR()`, etc.), we do not assume that all agents in the model are connected with each other. Thus, a network of agents must be built using the `agents_smallworld()` function before running the model where:\
- n = number of agents\
- k = number of ties in the small world network\
- d = whether the graph is directed or not\
- p = probability of rewiring
:::
## Important Statistics
```{r get-hist-total}
head(get_hist_total(model_sirconn))
```
An important statistic in epidemiological models is the reproductive number.
```{r repnum}
repnum <- get_reproductive_number(model_sirconn)
head(repnum)
```
**`epiworldR`** has a method to automatically plot the reproductive number. Thisfunction takes the average of values in the above table for each date and repeats until all date have been accounted for. For example, on average, individuals who acquired the virus on the 10th day transmit the virus to roughly 1.7 other individuals.
```{r}
x <- plot(repnum, type="b")
```
## Exercise 1
Create a SEIR model using the `ModelSEIR()` function (not `ModelSEIRCONN()`) to simulate a COVID-19 outbreak for 100 days in a population with:\
- prevalence = 0.01\
- transmission_rate = 0.9\
- recovery_rate = $\frac{1}{4}$\
- incubation_days = 4
Then plot the model parameters to analyze changes in counts over time. When running the model, set seed = 1912.
To accomplish this for a SEIR model, you will need to add the model to a smallworld population using the `agents_smallworld()` function after initializing the model. From there, run the model and visualize. Assume:\
- n = 10000\
- k = 5\
- d = FALSE\
- p = .01
After how many days does the number of infections peak in this simulation? How many infections occur at the peak?
```{r}
# Your solution here
```
## Exercise 2
Plot the reproductive number of the COVID-19 simulated SEIR model over 100 days.
```{r}
# Your solution here
```