In epiffiter
, opposite to the fit_
functions (estimate parameters from fitting models to the data), the sim_
family of functions allows to produce the DPC data given a set of parameters for a specific model. Currently, the same four population dynamic models that are fitted to the data can be simulated.
The functions use the ode()
function of the devolve
package (Soetaert,Petzoldt & Setzer 2010) to solve the differential equation form of the e epidemiological models.
First, we need to load the packages we'll need for this tutorial.
library(epifitter) library(magrittr) library(ggplot2) library(cowplot)
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 replicatesWhen 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 replicatesLet'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)
A 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 noiseUse 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
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)
mono_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
logist_model <- sim_logistic( N = 100, y0 = 0.01, dt = 5, r = 0.1, alpha = 0.2, n = 7 ) head(logist_model)
logist_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
gomp_model <- sim_gompertz( N = 100, y0 = 0.01, dt = 5, r = 0.07, alpha = 0.2, n = 7 ) head(gomp_model)
gomp_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
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)
Karline Soetaert, Thomas Petzoldt, R. Woodrow Setzer (2010). Solving Differential Equations in R: Package deSolve. Journal of Statistical Software, 33(9), 1--25. DOI: 10.18637/jss.v033.i09
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.