Getting Started in R¶

Prerequisites¶

Open R or RStudio.

Make sure you have installed MetaWards according to the R installation instructions. Check this is the case by typing;

> metawards::py_metawards_available()
[1] TRUE


If you don’t see TRUE returned, then double-check your installation.

Now make sure that you have installed the tidyverse, e.g. via;

> install.packages("tidyverse")
> library(tidyverse)


Then, finally, load the metawards R module;

> library(metawards)


Creating the disease¶

You should now be in R (or RStudio) and have imported metawards. To run a simulation you need to define the Disease that you want to model. MetaWards implements a SEIR-style model, but you have complete control to define as many (or few) stages as you wish.

First, we will create a disease, which we will call lurgy, that will consist of four stages: S, E, I and R. To do this, let’s create the disease;

> lurgy <- metawards$Disease(name="lurgy")  Next, we will add each stage. You don’t define the “S” stage, as the model starts with a set of susceptible individuals by default. Instead, we need to add in the E, I and R stages. First, lets add the latent (“E”) stage. Latent individuals are not infectious, and so we will set beta (the infectivity parameter) to 0.0. Individuals will progress quickly through this stage, so we will set progress to 0.5, meaning that 50% of individuals move to the next stage each day. > lurgy$add("E", beta=0.0, progress=0.5)


Next we will add the infectious (“I”) stage. This will have a high beta value (0.8), but a lower progress (0.25) as we will model this as a disease with a long symptomatic period.

> lurgy$add("I", beta=0.8, progress=0.25)  Finally, we need to add the recovered (“R”) stage. We don’t need to set the beta or progress values, as MetaWards will automatically recognise this as the recovered state, and will set beta to 0 and progress to 0 automatically. > lurgy$add("R")


You can should print this disease to the screen to confirm that everything has been correctly set.

> print(lurgy)

* Disease: lurgy
* stage: ['E', 'I', 'R']
* mapping: ['E', 'I', 'R']
* beta: [0, 0.8, 0]
* progress: [0.5, 0.25, 0]
* too_ill_to_move: [0, 0, 0]
* start_symptom: 2


Note

You can save this disease to a file using lurgy$to_json("lurgy.json.bz2"), and then load it back using lurgy = metawards$Disease$from_json("lurgy.json.bz2") Creating the wards (network)¶ Next, you need to define the wards (network) that will contain the individuals who will experience the model outbreak. We will first start with a single ward, called home. > home <- metawards$Ward(name="home")


MetaWards works by assigning individuals as either workers or players. The difference is that workers make fixed (predictable) movements between different wards each day, while players make random movements. Since we have just a single ward, we will start by populating it with 10,000 players.

> home$set_num_players(10000) > print(home) Ward( info=home, num_workers=0, num_players=10000 )  Note You can save this Ward to a file using home$to_json("home.json.bz2"), and then load it back using home = metawards$Ward$from_json("home.json.bz2")

Running the model¶

Now we have a disease and a network, we can now model an outbreak. To do this, we will use the metawards.run() function.

> results <- metawards$run(model=home, disease=lurgy)  This will print a lot to the screen. The key lines are these; ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ Day 0 ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ S: 10000 E: 0 I: 0 R: 0 IW: 0 POPULATION: 10000 ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ Day 1 ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ S: 10000 E: 0 I: 0 R: 0 IW: 0 POPULATION: 10000 Number of infections: 0 ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ Day 2 ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ S: 10000 E: 0 I: 0 R: 0 IW: 0 POPULATION: 10000 Number of infections: 0 ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ Day 3 ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ S: 10000 E: 0 I: 0 R: 0 IW: 0 POPULATION: 10000 Number of infections: 0 ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ Day 4 ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ S: 10000 E: 0 I: 0 R: 0 IW: 0 POPULATION: 10000 Number of infections: 0 ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ Day 5 ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ S: 10000 E: 0 I: 0 R: 0 IW: 0 POPULATION: 10000 Number of infections: 0 Ending on day 5  This shows the number of people in the different stages of the outbreak. In this case, there was no infection seeded, and so the number of infections remained zero. Seeding the outbreak¶ We need to seed the outbreak with some additional seeds. We do this using the additional option. This can be very powerful (e.g. adding seeds at different days, different wards etc.), but at its simplest, it is just the number of initial infections on the first day in the first ward. We will start with 100 initial infections; > results <- metawards$run(model=home, disease=lurgy, additional=100)


Now you get a lot more output, e.g. for me the outbreak runs for 75 days.

━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ Day 70 ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
S: 423  E: 0  I: 1  R: 9576  IW: 0  POPULATION: 10000
Number of infections: 1

━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ Day 71 ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
S: 423  E: 0  I: 1  R: 9576  IW: 0  POPULATION: 10000
Number of infections: 1

━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ Day 72 ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
S: 423  E: 0  I: 1  R: 9576  IW: 0  POPULATION: 10000
Number of infections: 1

━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ Day 73 ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
S: 423  E: 0  I: 1  R: 9576  IW: 0  POPULATION: 10000
Number of infections: 1

━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ Day 74 ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
S: 423  E: 0  I: 1  R: 9576  IW: 0  POPULATION: 10000
Number of infections: 1

━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ Day 75 ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
S: 423  E: 0  I: 0  R: 9577  IW: 0  POPULATION: 10000
Number of infections: 0
Ending on day 75


Visualising the results¶

The output results contains the filename of a csv file that contains the S, E, I and R data (amongst other things). You can load and plot this using standard R commands, e.g.

> results <- read.csv(results)
> print(results)
fingerprint repeat. day       date     S    E    I    R IW SCALE_UV
1       REPEAT       1   0 2020-07-20 10000    0    0    0  0        1
2       REPEAT       1   1 2020-07-21  9900   57   43    0  1        1
3       REPEAT       1   2 2020-07-22  9859   66   66    9  1        1
4       REPEAT       1   3 2020-07-23  9807   86   82   25  1        1
5       REPEAT       1   4 2020-07-24  9747  101  112   40  1        1
6       REPEAT       1   5 2020-07-25  9654  140  130   76  1        1
7       REPEAT       1   6 2020-07-26  9548  183  165  104  1        1
8       REPEAT       1   7 2020-07-27  9433  215  203  149  1        1
9       REPEAT       1   8 2020-07-28  9280  252  269  199  1        1
10      REPEAT       1   9 2020-07-29  9082  318  341  259  1        1
...


To visualise the data we need to tidy it up so that we can group by S, E, I and R.

> results <- results %>%
pivot_longer(c("S", "E", "I", "R"),
names_to = "stage", values_to = "count")
> print(results)
# A tibble: 304 x 8
fingerprint repeat.   day date          IW SCALE_UV stage count
<chr>         <int> <int> <chr>      <int>    <dbl> <chr> <int>
1 REPEAT            1     0 2020-07-20     0        1 S     10000
2 REPEAT            1     0 2020-07-20     0        1 E         0
3 REPEAT            1     0 2020-07-20     0        1 I         0
4 REPEAT            1     0 2020-07-20     0        1 R         0
5 REPEAT            1     1 2020-07-21     1        1 S      9900
6 REPEAT            1     1 2020-07-21     1        1 E        57
7 REPEAT            1     1 2020-07-21     1        1 I        43
8 REPEAT            1     1 2020-07-21     1        1 R         0
9 REPEAT            1     2 2020-07-22     1        1 S      9859
10 REPEAT            1     2 2020-07-22     1        1 E        66
# … with 294 more rows


You can graph S, E, I and R against day using;

> ggplot(data = results,
mapping = aes(x=day, y=count, color=stage)) + geom_line()


The result should look something like this;

Complete code¶

The complete R code for this part of the getting started guide is re-copied below;

# Load the dependencies / libraries
library(tidyverse)
library(metawards)

# Create the disease
lurgy <- metawards$Disease(name="lurgy") lurgy$add("E", beta=0.0, progress=0.25)
lurgy$add("I", beta=0.8, progress=0.25) lurgy$add("R")

# Create the model network
home <- metawards$Ward(name="home") home$set_num_players(10000)

# Run the model