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)
devtools::load_all()
```
## 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}
sim <- seiqhrf(init, control, param)
sim
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
```
### Extract summary of the simulation
```{r, extract summary}
res <- summary(sim)
names(res)
```
### 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(sim, "times")
```
### Plot baseline sirplus results
To visualise your sirplus model, you can plot the change in prevalence (i.e. people)
over time in each compartment by changing parameter `method` to `models`.
```{r viz prevalence}
plot(sim,
start_date = lubridate::ymd("2020-01-01"),
comp_remove = c('s.num', 'r.num'),
plot_title = 'Baseline Model')
```