Hands On
Packages
First, we need to load the packages we’ll need for this tutorial.
library(epifitter)
library(magrittr)
library(ggplot2)
library(cowplot)The basics of the simulation
The sim_ functions, regardless of the model, require the same set of six arguments. By default, at least two arguments are required (the others have default values)
r: apparent infection raten: number of replicates
When n is greater than one, replicated epidemics (e.g. replicated treatments) are produced and a level of noise (experimental error) should be set in the alpha argument. These two arguments combined set will generate random_y values, which will vary randomly across the defined number of replicates.
The other arguments are:
N: epidemic duration in time unitsdt: time (fixed) in units between two assessmentsy0: initial inoculumalpha: noise parameters for the replicates
Exponential
Let’s simulate a curve resembling the exponential growth.
exp_model <- sim_exponential(
N = 100, # total time units
y0 = 0.01, # initial inoculum
dt = 10, # interval between assessments in time units
r = 0.045, # apparent infection rate
alpha = 0.2,# level of noise
n = 7 # number of replicates
)
head(exp_model)## replicates time y random_y
## 1 1 0 0.01000000 0.01412441
## 2 1 10 0.01568425 0.01729569
## 3 1 20 0.02459905 0.01934236
## 4 1 30 0.03858028 0.04199333
## 5 1 40 0.06050749 0.04437351
## 6 1 50 0.09489670 0.09857329A data.frame object is produced with four columns:
replicates: the curve with the respective ID numbertime: the assessment timey: the simulated proportion of disease intensityrandom_y: randomly simulated proportion disease intensity based on the noise
Use the ggplot2 package to build impressive graphics!
exp_plot = exp_model %>%
ggplot(aes(time, y)) +
geom_jitter(aes(time, random_y), size = 3,color = "gray", width = .1) +
geom_line(size = 1) +
theme_minimal_hgrid() +
ylim(0,1)+
labs(
title = "Exponential",
y = "Disease intensity",
x = "Time"
)
exp_plot
Monomolecular
The logic is exactly the same here.
mono_model <- sim_monomolecular(
N = 100,
y0 = 0.01,
dt = 5,
r = 0.05,
alpha = 0.2,
n = 7
)
head(mono_model)## replicates time y random_y
## 1 1 0 0.0100000 0.0100000
## 2 1 5 0.2289861 0.2767016
## 3 1 10 0.3995322 0.4386150
## 4 1 15 0.5323535 0.5588012
## 5 1 20 0.6357949 0.6653795
## 6 1 25 0.7163551 0.7211241mono_plot = mono_model %>%
ggplot(aes(time, y)) +
geom_jitter(aes(time, random_y), size = 3, color = "gray", width = .1) +
geom_line(size = 1) +
theme_minimal_hgrid() +
labs(
title = "Monomolecular",
y = "Disease intensity",
x = "Time"
)
mono_plot
The Logistic model
logist_model <- sim_logistic(
N = 100,
y0 = 0.01,
dt = 5,
r = 0.1,
alpha = 0.2,
n = 7
)
head(logist_model)## replicates time y random_y
## 1 1 0 0.01000000 0.01000000
## 2 1 5 0.01638216 0.01517832
## 3 1 10 0.02672677 0.02864972
## 4 1 15 0.04331509 0.04849382
## 5 1 20 0.06946352 0.04272364
## 6 1 25 0.10958806 0.09518275logist_plot = logist_model %>%
ggplot(aes(time, y)) +
geom_jitter(aes(time, random_y), size = 3,color = "gray", width = .1) +
geom_line(size = 1) +
theme_minimal_hgrid() +
labs(
title = "Logistic",
y = "Disease intensity",
x = "Time"
)
logist_plot
Gompertz
gomp_model <- sim_gompertz(
N = 100,
y0 = 0.01,
dt = 5,
r = 0.07,
alpha = 0.2,
n = 7
)
head(gomp_model)## replicates time y random_y
## 1 1 0 0.01000000 0.01236849
## 2 1 5 0.03896283 0.05242822
## 3 1 10 0.10158896 0.09898759
## 4 1 15 0.19958740 0.15206340
## 5 1 20 0.32122825 0.31599946
## 6 1 25 0.44922018 0.49399206gomp_plot = gomp_model %>%
ggplot(aes(time, y)) +
geom_jitter(aes(time, random_y), size = 3,color = "gray", width = .1) +
geom_line(size = 1) +
theme_minimal_hgrid() +
labs(
title = "Gompertz",
y = "Disease intensity",
x = "Time"
)
gomp_plot
Combo
Use the function plot_grid() from the cowplot package to gather all plots into a grid
plot_grid(exp_plot,
mono_plot,
logist_plot,
gomp_plot)