Newer
Older
---
title: "Introduction to sirplus (V2)"
date: "Last updated: 12 April 2020"
vignette: >
%\VignetteIndexEntry{Introduction to sirplus (V2)}
%\VignetteEngine{knitr::rmarkdown}
\usepackage[utf8]{inputenc}
---
```{r setup, include=FALSE}
version_date <- lubridate::ymd("2020-04-12")
knitr::opts_chunk$set(echo = TRUE, cache=FALSE, eval=TRUE, collapse = TRUE,
tidy.opts=list(width.cutoff=60), tidy=FALSE)
```
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)
```
## Simulate and inspect a baseline sirplus model
### Set parameters
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, set parameters}
s.num <- 2000 # number susceptible
i.num <- 15 # number infected
q.num <- 5 # number in self-isolation
h.num <- 1 # number in the hospital
nsteps <- 90 # number of steps (e.g. days) to simulate
control <- control_seiqhrf(nsteps = nsteps)
init <- init_seiqhrf(s.num = s.num, i.num = i.num, q.num = q.num, h.num = h.num)
print(init)
print(control)
print(param)
```
### Simulate baseline
This will produce an seiqhrf object.
```{r, run simulation}
```
### Extract summary of the simulation
```{r, extract summary}
```
### 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 `plot()` functions, and set parameter `method` as `times`
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 baseline sirplus results
To visualise your sirplus model, you can plot the change in prevalence (i.e. people)
```{r viz prevalence}
start_date = lubridate::ymd("2020-01-01"),
comp_remove = c('s.num', 'r.num'),
plot_title = 'Baseline Model')
```
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
## 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 = nsteps, vals = vals, timing = timing)
control <- control_seiqhrf(nsteps = nsteps)
param <- param_seiqhrf(act.rate.e = act_rate, act.rate.i = act_rate * 0.5)
init <- init_seiqhrf(s.num = s.num, i.num = i.num, q.num = q.num, h.num = h.num)
sim_exp <- seiqhrf(init, control, param)
# Experiment #2
vals <- c(10, 7, 7, 10)
timing <- c(7, 14, 21, 28)
act_rate_relax <- vary_param(nstep = nsteps, vals = vals, timing = timing)
param <- param_seiqhrf(act.rate.e = act_rate_relax, act.rate.i = act_rate_relax * 0.5)
sim_exp_relax <- seiqhrf(init, control, param)
# Compare experiments 1 and 2 to the baseline simulation
plot(list("Baseline" = baseline_sim, "Closures" = sim_exp, "Closures (2 mo)" = sim_exp_relax),
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
Use the `plot()` function with parameter `method = "weekly_local"` to extract and plot the number of weekly expected patients in a hospital. If `return_df` is `TRUE`, the function returns a dataframe that generates the plot, otherwise, it returns a ggplot object.
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_lim`: 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}
plot(baseline_sim, method = "weekly_local", time_lim = 50, start_date = lubridate::ymd("2020-01-01"), return_df = TRUE)
```