---
title: "Getting started with epiworldR"
author:
  - Derek Meyer
  - George Vega Yon
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Getting started with epiworldR}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>", out.width = "80%", fig.width = 7, fig.height = 5,
  fig.align = "center"
)
```



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-virus 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.)

# Example 1: Simulating an SIR model

## Setup and running the model

This example implements a social network with parameters listed within the
`ModelSIRCONN` function. The virus name is specified (COVID-19), 50000 agents
are initialized, the virus prevalence of 0.001 is declared, each agent will
contact two others (contact_rate), the transmission rate for
any given agent is 0.3, and the recovery rate is set to $\frac{1}{3}$.
To create this model on epiworldR, simply use the `ModelSIRCONN` function.
From here, the example will take you through the basic features of epiworldR.

```{r sirconn-setup}
library(epiworldR)
model_sir <- ModelSIRCONN(
  name              = "COVID-19",
  n                 = 50000,
  prevalence        = 0.0001,
  contact_rate      = 2,
  transmission_rate = 0.5,
  recovery_rate     = 1 / 3
)

# Printing the model
model_sir
```

Printing the model shows us some information. Nevertheless, we can extract detailed information using the summary method.

```{r summary-method}
summary(model_sir)
```

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 viruses, amount of time the last replicate took to run (last run elapsed t), and rewiring
status (on or off). The model also includes a list of global actions (interventions) that are called during the model run. Next, you will see a list of the
viruses used in the model. In this case, COVID-19 was the only virus used.
Note that epiworldR can include more than one virus in a
model. Tool(s) lists agents' tools to fight the virus. Examples of
this may include masking, vaccines, social distancing, etc. In this model, no
tools are specified. Lastly, model parameters are listed.


To execute the model, use the run function with the SIR model object, the 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_sir, ndays = 50, seed = 1912)
summary(model_sir)
```

```{r getting-totals, echo=FALSE}
initials <- get_hist_total(model_sir)[1:3, ]$counts |> prettyNum(big.mark = ",")
finals   <- get_today_total(model_sir) |> prettyNum(big.mark = ",")
tmat     <- get_transition_probability(model_sir)
tmat     <- round(tmat, digits = 2)
```

There are two additional sections included in the summary after running the model. First, we see the distribution of the population at time 50. This section describes the flow of agents from each state (SIR) after 50 days. In the
example, you'll see the number of agents in the susceptible state decreased
from `r initials[1]` to `r finals[1]`, the number of agents in the infected state increased from `r initials[2]` to `r finals[2]`, and recovered agents increased to `r finals[3]` after 50 days.
The counts for these states will 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. For example,
in the susceptible row, each agent has a `r tmat[1]` probability of remaining in the susceptible state with a `r tmat[1,2]` probability of moving from the susceptible state to the infected state. Notice there is no chance of skipping states. In other words, an agent can't jump from the
susceptible state to the recovered state; that agent must pass through the
infected state to progress to the recovered state. The same logic
applies to moving backward; an agent cannot become susceptible again after
infection.

## Plot

## Extracting information

After running the epiworldR model, below is a list of all the functions that can
be called using the epiworld model object.

```{r showing-methods}
methods(class = "epiworld_model")
```

To demonstrate, start with the basic plot and get_hist_total functions.

```{r}
plot(model_sir)
```

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.

```{r get-hist-total}
head(get_hist_total(model_sir))
```

An essential statistic in epidemiological models is the reproductive number:

```{r repnum}
repnum <- get_reproductive_number(model_sir)
head(repnum)
```

epiworldR has a method to plot the reproductive number automatically.
The function takes the average of values in the above table for each date and
repeats until all data have been accounted for.

```{r}
x <- plot(repnum, type = "b")
subset(x, date == 10) # Reproductive number on day 10
```

Another typical piece of information is the daily incidence. This is the number of new
cases per day. In epiworldR, we can get the incidence by looking at the daily
transitions between states. Although the function `get_hist_transition_matrix` provides the desired data, the function `plot_incidence` is a nice wrapper for visualizing the data:

```{r}
plot_incidence(model_sir)
```

## Adding more viruses/viruses

epiworldR supports multi-virus models. The below code gives instructions on
how to implement this. Using the `virus` function, give a name to the new
virus/virus with its corresponding probability of infecting any given agent.
In this example, `prob_infecting` is set to 1.0, making it highly contagious.
To officially add this new virus to the model, use the `add_virus`
function by calling the original epiworldR model object, the new virus, and
the new virus' prevalence (which is set to 0.01 in this example).

```{r design-and-add}
# Building the virus
flu <- virus(
  name = "Flu", prob_infecting = .3,
  prevalence = .0001, as_proportion = TRUE
)

# Adding the virus to the model
add_virus(model_sir, flu)
```

After running the updated model with the new virus included for 50 days, the
output below describes the simulation. To confirm that the flu is included,
notice the presence of "Flu" in the Virus(es) section of the output. All other
output is interpretable as specified in previous sections.

```{r}
run(model_sir, ndays = 50, seed = 1912)
model_sir
```

Plotting the previous model (including the flu) yields the following. Notice
the presence of two reproductive numbers plotted over time. Variant 0 refers
to COVID-19, and virus 1 refers to the flu.

```{r, fig.height=10}
repnum2 <- get_reproductive_number(model_sir)

op <- par(mfrow = c(2, 1))
plot(model_sir)
plot(repnum2, type = "b")
par(op)
```

## Tools

Now, the implementation of tools to combat any viruses and viruses in the model
will be demonstrated. First, for the sake of simplicity, remove the flu virus
from the SIR model object (keep in mind the index for the flu virus in the
model object is 1). Next, provide parameters for the new tool using the `tool`
function. These parameters include the name of the tool, any reduction in
probabilities for the SIR model parameters, and increased probability of
recovery option. In order to add the tool to the SIR model, use the `add_tool`
function with the SIR model object, new tool, and prevalence of the tool.
In this example, assume that 85% of the population will have received the
vaccination.

```{r}
# Removing the flu virus from the model
rm_virus(model_sir, 1)

vaccine <- tool(
  name = "Vaccine",
  prevalence = 0.5,
  as_proportion = TRUE,
  susceptibility_reduction = .9,
  transmission_reduction = .5,
  recovery_enhancer = .5,
  death_reduction = .9
)

add_tool(model_sir, vaccine)
run(model_sir, ndays = 50, seed = 1231)
```

```{r curves-including-vaccine, fig.height=10}
repnum3 <- get_reproductive_number(model_sir)

op <- par(mfrow = c(2, 1))
plot_incidence(model_sir)
plot(repnum3, type = "b")
par(op)
```