---
title: "Introduction to sirplus"
date: "Last updated: 26 March 2020"
vignette: >
  %\VignetteIndexEntry{Introduction to sirplus}
  %\VignetteEngine{knitr::rmarkdown}
  \usepackage[utf8]{inputenc}
---

```{r setup, include=FALSE, eval=TRUE}
version_date <- lubridate::ymd("2020-03-26")
knitr::opts_chunk$set(echo = TRUE, cache=FALSE, eval=TRUE, collapse = TRUE,
                      tidy.opts=list(width.cutoff=60), tidy=FALSE)
```

![](sirplus-logo-small.png)

The sirplus package makes it easy to generate stochastic individual compartment 
models (ICMs) to simulate contagious disease spread using compartments not 
available in standard SIR packages. This method and most of the code was 
originally written by Tim Churches (see his [blog post](https://timchurches.github.io/blog/posts/2020-03-18-modelling-the-effects-of-public-health-interventions-on-covid-19-transmission-part-2/)). The sirplus package was developed by 
the Bioinformatics & Cellular Genomics team at St. Vincent's Institute of 
Medical Research in order to help St. Vincents' Hospital model the COVID-19
pandemic.

The compartments available in this package include: 

| Compartment | Functional definition                                                             |
|-------------|-----------------------------------------------------------------------------------|
| S           | Susceptible individuals                                                           |
| E           | Exposed **and** infected, not yet symptomatic but potentially infectious          |
| I           | Infected, symptomatic **and** infectious                                          |
| Q           | Infectious, but (self-)isolated                                                   |
| H           | Requiring hospitalisation (would normally be hospitalised if capacity available)  |
| R           | Recovered, immune from further infection                                          |
| F           | Case fatality (death due to infection, not other causes)                          |

Again, for more information about these compartments and about the parameters 
that are associated with them see Tim Churches's [blog post](https://timchurches.github.io/blog/posts/2020-03-18-modelling-the-effects-of-public-health-interventions-on-covid-19-transmission-part-2/).

```{r Load package}
#library(sirplus)
devtools::load_all(".")
```

## Simulate and inspect a baseline sirplus model

### Simulate

Here we will simulate the epidemiological data for a made-up population with 1000 
susceptible individuals (S), 50 that are infected but not in the hospital or in 
self-quarantine (I; maybe people that are infected/symptomatic but not tested/
aware), 10 confirmed cases that have self-isolated (Q), and 1 confirmed case that
has been hospitalized (H). We call this the baseline model because it uses 
default parameters for disease spread (i.e. no additional interventions).

```{r simulate baselines}  
s.num <- 2000  # number susceptible
i.num <- 15  # number infected 
q.num <- 5  # number in self-isolation
h.num <- 1  # number in the hospital
nstep <- 90 # number of steps (e.g. days) to simulate

baseline_sim <- simulate_seiqhrf(nstep = nstep, s.num = s.num, i.num = i.num,
                                 q.num = q.num, h.num = h.num)

head(baseline_sim$df)
```

### Inspect baseline transition distributions

The sirplus model controls transitions between compartments, i.e. a change in state for an individual (e.g. going from self-isolation to hospital), using a variety of 
transition parameters. You can use the `get_times()` and `plot_times()` 
functions to examine the distributions of timings for various transitions based
on these parameters. In the case of a disease with observed data available, these plots can be used to sanity check parameter settings.

```{r baseline sims, fig.height=8}
plot_times(baseline_sim$sim)
```


### Plot baseline sirplus results

To visualise your sirplus model, you can plot the change in prevalence (i.e. people) 
over time in each compartment.

```{r viz prevalence}
plot_models(sims = baseline_sim,
             start_date = lubridate::ymd("2020-01-01"),
             comp_remove = c('s.num', 'r.num'),
             plot_title = 'Baseline Model')
```

## Run an experiment

With the sirplus package you can also set up experiments. We will set up two
experiments here: 

- Experiment #1: One week after the beginning of the epidemic, schools and non-essential businesses are closed to encourage social distancing. This causes the act.rate to gradually drop from 10 to 6 over the course of the next week. In this experiment, we imagine these policies are never lifted, so act.rate remains at 6 for the duration of the simulation. 
- Experiment #2: Again, one week after the beginning of the epidemic, social distancing policies are put into place resulting in act.rate dropping from 10 to 6 over the next week. But after two weeks these policies are lifted and the act.rate returns to normal within the next week. 

```{r Experiment example with act.rate}
# Experiment #1
vals <- c(10, 7)
timing <- c(7, 14)
act_rate <- vary_param(nstep = nstep, vals = vals, timing = timing)

sim_exp <- simulate_seiqhrf(nstep = nstep, s.num = s.num, i.num = i.num,
                            q.num = q.num, h.num = h.num, act.rate.e = act_rate, 
                            act.rate.i = act_rate * 0.5)

# Experiment #2
vals <- c(10, 7, 7, 10)
timing <- c(7, 14, 21, 28)
act_rate_relax <- vary_param(nstep, vals = vals, timing = timing)

sim_exp_relax <- simulate_seiqhrf(nstep = nstep, s.num = s.num, i.num = i.num,
                            q.num = q.num, h.num = h.num, 
                            act.rate.e = act_rate_relax, 
                            act.rate.i = 0.5 * act_rate_relax)

# Compare experiments 1 and 2 to the baseline simulation
plot_models(sims = c(baseline_sim, sim_exp, sim_exp_relax),
            sim_id = c('Baseline', 'Closures', 'Closures (2 mo)'),
            start_date = lubridate::ymd("2020-01-01"),
            comp_remove = c('s.num', 'r.num'),
            plot_title = 'Closures Experiment')

```

From these results we see that this policy would likey reduce the peak number
of infections from 500 to 300 and would "flatten the curve" for infections and 
thus hospitalizations (the peak in number of infected people is lower, but the number of people infected declines more slowly after the peak of the epidemic than in the baseline model).


## Visualize sirplus models to aid hospitals

Use the `get_weekly_local()` function to extract and plot the number of weekly expected 
patients in a hospital. This function takes the following input:

  - `sims`: Single or list of outputs from simulate.seiqhrf().
  - `market.share`: Percent of cases in your model population (s.num) that are
  anticipated at the hospital of interest. Default = 4% (0.04).
  - `icu_percent`: Percent of hospitalised cases are are likely to need treatment in an intensive care unit (ICU). Default = 10% (0.1).
  - `start_date`: Epidemic start date. Default is `NA`.
  - `time_limit`: Number of days (nsteps) to include. Default = 90.
  - `total_population`: True population size. This parameter is only needed if
  the simulation size (s.num) was smaller than the true population size (i.e. 
  scaled down) to reduce computational cost.
  
```{r hospital visualization}
weekly_local <- get_weekly_local(baseline_sim$sim, time_lim = 50, start_date = lubridate::ymd("2020-01-01"))

head(weekly_local$result)
weekly_local$plot
```
