Examples

This R package includes several popular epidemiological models, including SIS (wiki), SIR (wiki), and SEIR (wiki) using either a fully connected graph (similar to a compartmental model) or a user-defined network. Here are some examples:

SIR Model Using a Random Graph

This Susceptible-Infected-Recovered model features a population of 100,000 agents simulated in a small-world network. Each agent is connected to ten other agents. One percent of the population has the virus, with a 70% chance of transmission. Infected individuals recover at a 0.3 rate:

Code

library(epiworldR)

## Creating a SIR model
sir <- ModelSIR(
  name           = "COVID-19",
  prevalence     = .01,
  transmission_rate  = .7,
  recovery_rate =  .3
  ) |>
  # Adding a Small world population 
  agents_smallworld(n = 100000, k = 10, d = FALSE, p = .01) |>
  # Running the model for 50 days
  run(ndays = 50, seed = 1912)

_________________________________________________________________________
|Running the model...
|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| done.
| done.

Code

sir

________________________________________________________________________________
Susceptible-Infected-Recovered (SIR)
It features 100000 agents, 1 virus(es), and 0 tool(s).
The model has 3 states.
The final distribution is: 822 Susceptible, 415 Infected, and 98763 Recovered.

Visualizing the outputs

Code

plot(sir)

SEIR Model With a Fully Connected Graph

Code

model_seirconn <- ModelSEIRCONN(
  name                = "COVID-19",
  prevalence          = 0.01, 
  n                   = 10000,
  contact_rate        = 4, 
  incubation_days     = 7, 
  transmission_rate   = 0.6,
  recovery_rate       = 0.5
)

set.seed(132)
run(model_seirconn, ndays = 100)

_________________________________________________________________________
Running the model...
||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| done.
 done.

Code

model_seirconn

________________________________________________________________________________
Susceptible-Exposed-Infected-Removed (SEIR) (connected)
It features 10000 agents, 1 virus(es), and 0 tool(s).
The model has 4 states.
The final distribution is: 91 Susceptible, 0 Exposed, 0 Infected, and 9909 Recovered.

Computing some key statistics: plotting and reproductive number (wiki)

Code

plot(model_seirconn)

Code

repnum <- get_reproductive_number(model_seirconn)
plot(repnum, type = "b")

SIR Logit

Code

set.seed(2223)
n <- 100000

X <- cbind(
  Intercept = 1,
  Female    = sample.int(2, n, replace = TRUE) - 1
  )

coef_infect  <- c(.1, -2, 2)
coef_recover <- rnorm(2)

model_logit <- ModelSIRLogit(
  "covid2",
  data = X,
  coefs_infect      = coef_infect,
  coefs_recover     = coef_recover, 
  coef_infect_cols  = 1L:ncol(X),
  coef_recover_cols = 1L:ncol(X), 
  prob_infection = .8,
  recovery_rate = .3,
  prevalence = .01
)

agents_smallworld(model_logit, n, 8, FALSE, .01)

run(model_logit, 50)

_________________________________________________________________________
|Running the model...
|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| done.
| done.

Code

plot(model_logit)

Code

## Females are supposed to be more likely to become infected
rn <- get_reproductive_number(model_logit)

(table(
  X[, "Female"],
  (1:n %in% rn$source)
) |> prop.table())[,2]

      0       1 
0.12984 0.14201 

Code

## Looking into the agents
get_agents(model_logit)

Agents from the model "Susceptible-Infected-Removed (SIR) (logit)":
Agent: 0, state: Recovered (2), Nvirus: 0, NTools: 0, NNeigh: 8
Agent: 1, state: Recovered (2), Nvirus: 0, NTools: 0, NNeigh: 8
Agent: 2, state: Recovered (2), Nvirus: 0, NTools: 0, NNeigh: 8
Agent: 3, state: Recovered (2), Nvirus: 0, NTools: 0, NNeigh: 8
Agent: 4, state: Recovered (2), Nvirus: 0, NTools: 0, NNeigh: 8
Agent: 5, state: Recovered (2), Nvirus: 0, NTools: 0, NNeigh: 8
Agent: 6, state: Recovered (2), Nvirus: 0, NTools: 0, NNeigh: 8
Agent: 7, state: Recovered (2), Nvirus: 0, NTools: 0, NNeigh: 8
Agent: 8, state: Susceptible (0), Nvirus: 0, NTools: 0, NNeigh: 8
Agent: 9, state: Recovered (2), Nvirus: 0, NTools: 0, NNeigh: 8
... 99990 more agents ...

Transmission Network

Code

sir <- ModelSIR(
  name              = "COVID-19",
  prevalence        = .01,
  transmission_rate = .5,
  recovery_rate     = .5
  ) |>
    # Adding a Small world population 
    agents_smallworld(n = 500, k = 10, d = FALSE, p = .01) |>
    # Running the model for 50 days
    run(ndays = 50, seed = 1912)

_________________________________________________________________________
|Running the model...
|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| done.
| done.

Code

## Transmission network
net <- get_transmissions(sir)

## Plotting
library(netplot)

Loading required package: grid

Code

library(igraph)

Attaching package: 'igraph'

The following object is masked from 'package:netplot':

    ego

The following objects are masked from 'package:stats':

    decompose, spectrum

The following object is masked from 'package:base':

    union

Code

x <- graph_from_edgelist(as.matrix(net[,2:3]) + 1)

nplot(x, edge.curvature = 0, edge.color = "gray", skip.vertex=TRUE)

Multiple Simulations

Code

model_sir <- ModelSIRCONN(
  name              = "COVID-19",
  prevalence        = 0.01,
  n                 = 1000,
  contact_rate      = 2,
  transmission_rate = 0.9,
  recovery_rate     = 0.1
  )

## Generating a saver
saver <- make_saver("total_hist", "reproductive")

## Running and printing
run_multiple(model_sir, ndays = 100, nsims = 50, saver = saver, nthread = 2)

Starting multiple runs (50) using 2 thread(s)
_________________________________________________________________________
_________________________________________________________________________
||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| done.
 done.

Code

## Retrieving the results
ans <- run_multiple_get_results(model_sir)

head(ans$total_hist)

  sim_num date nviruses       state counts
1       1    0        1 Susceptible    990
2       1    0        1    Infected     10
3       1    0        1   Recovered      0
4       1    1        1 Susceptible    974
5       1    1        1    Infected     25
6       1    1        1   Recovered      1

Code

head(ans$reproductive)

  sim_num virus_id    virus source source_exposure_date rt
1       1        0 COVID-19    767                   11  0
2       1        0 COVID-19    835                   10  0
3       1        0 COVID-19    466                    9  0
4       1        0 COVID-19    612                    9  0
5       1        0 COVID-19    793                    9  0
6       1        0 COVID-19     20                    8  0

Code

plot(ans$reproductive)