library(magrittr)
library(reshape2)
library(ggplot2)
library(deSolve)
## Warning: package 'deSolve' was built under R version 3.1.2
PlotFun <- function(parms = c(beta = (0.5*1), r = 1/7, u = 1/(70*365.25), v = 0.75),
inits = c(S = 499, I = 1, R = 0),
dt = seq(0, 52, 0.01)) {
fun <- function(t, x, parms) {
with(as.list(c(parms, x)), {
N <- S + I + R
dS <- - beta * I * S/N + (1-v) * u * N - u * S
dI <- + beta * I * S/N - r * I - u * I
dR <- + r * I + v * u * N - u * R
c(dS, dI, dR) %>% list
})
}
## simulate
simulation <- as.data.frame(lsoda(inits, dt, fun, parms = parms))
## Plot
outPlot <- melt(data = simulation,
id.vars = c("time"),
variable.name = "variable",
value.name = "value") %>%
subset(x = ., TRUE) %>%
ggplot(data = ., mapping = aes(x = time, y = value, color = variable)) +
layer(geom = "line", size = 2) +
theme_bw() + theme(legend.key = element_blank())
list(simulation, outPlot)
}
PlotFun(dt = seq(0, 52, 0.01))[[2]]
PlotFun(dt = seq(0, 100000, 10))[[2]]
PlotFun <- function(parms = c(beta = (0.5*1), r = 1/7, u = 1/(70*365.25), v = 0.75, e = 0.1),
inits = c(S = 499, I = 1, R = 0, V = 0),
dt = seq(0, 52, 0.01)) {
fun <- function(t, x, parms) {
with(as.list(c(parms, x)), {
N <- S + I + R + V
dS <- - beta * I * S/N + (1-v) * u * N - u * S
dI <- + beta * I * S/N + e*beta * I * V/N - r * I - u * I
dR <- + r * I - u * R
dV <- - e*beta * I * V/N + v * u * N - u * V
c(dS, dI, dR, dV) %>% list
})
}
## simulate
simulation <- as.data.frame(lsoda(inits, dt, fun, parms = parms))
## Plot
outPlot <- melt(data = simulation,
id.vars = c("time"),
variable.name = "variable",
value.name = "value") %>%
subset(x = ., TRUE) %>%
ggplot(data = ., mapping = aes(x = time, y = value, color = variable)) +
layer(geom = "line", size = 2) +
theme_bw() + theme(legend.key = element_blank())
list(simulation, outPlot)
}
PlotFun(dt = seq(0, 52, 0.01))[[2]]
PlotFun(dt = seq(0, 100000, 10))[[2]]
