Nothing
## ---- include = FALSE---------------------------------------------------------
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>"
)
## ----setup--------------------------------------------------------------------
library(MicroMoB)
library(ggplot2)
library(data.table)
library(parallel)
## -----------------------------------------------------------------------------
# mosquito parameters
f <- 0.3
q <- 1
eip <- 14
lifespan <- 20
g <- 1/lifespan
p <- 1- g
# human parameters
b <- 0.55
c <- 0.15
r <- 1/200
S <- 1e3
I <- 300
N <- S + I
# transmission parameters
kappa <- (I/N)*c
# equilibrium solutions
Z <- (r*I*N) / (b*f*q*S)
Y <- Z / (p^eip)
M <- (Z*(g + (f*q*p*kappa))) / (f*q*p*kappa*(p^eip))
lambda <- g*M
## -----------------------------------------------------------------------------
patches <- 1
nstrata <- 1
tmax <- 365 * 2
theta <- diag(nstrata)
psi <- diag(patches)
## -----------------------------------------------------------------------------
mod <- make_MicroMoB(tmax = tmax, p = patches)
setup_humans_SIS(mod, stochastic = FALSE, theta = theta, H = N, X = I, b = b, c = c, r = r)
setup_aqua_trace(mod, stochastic = FALSE, lambda = lambda)
setup_mosquito_RM(mod, stochastic = FALSE, f = f, q = q, eip = eip, p = p, psi = psi, M = M, Y = Y, Z = Z)
## -----------------------------------------------------------------------------
mosy_out <- data.table::CJ(day = 1:tmax, state = c('M', 'Y', 'Z'), value = NaN, species = "mosquito")
mosy_out <- mosy_out[c('M', 'Y', 'Z'), on="state"]
data.table::setkey(mosy_out, day)
human_out <- data.table::CJ(day = 1:tmax, state = c('S', 'I'), value = NaN, species = "human")
human_out <- human_out[c('S', 'I'), on="state"]
data.table::setkey(human_out, day)
while (get_tnow(mod) <= tmax) {
compute_bloodmeal_simple(model = mod)
step_aqua(model = mod)
step_mosquitoes(model = mod)
step_humans(model = mod)
mosy_out[day == get_tnow(mod) & state == 'M', value := mod$mosquito$M]
mosy_out[day == get_tnow(mod) & state == 'Y', value := mod$mosquito$Y]
mosy_out[day == get_tnow(mod) & state == 'Z', value := mod$mosquito$Z]
human_out[day == get_tnow(mod) & state == 'S', value := mod$human$H - mod$human$X]
human_out[day == get_tnow(mod) & state == 'I', value := mod$human$X]
mod$global$tnow <- mod$global$tnow + 1L
}
det_out <- rbind(mosy_out, human_out)
## -----------------------------------------------------------------------------
sto_out <- mclapply(X = 1:10, FUN = function(runid) {
mod <- make_MicroMoB(tmax = tmax, p = patches)
setup_humans_SIS(mod, stochastic = TRUE, theta = theta, H = N, X = I, b = b, c = c, r = r)
setup_aqua_trace(mod, stochastic = TRUE, lambda = lambda)
setup_mosquito_RM(mod, stochastic = TRUE, f = f, q = q, eip = eip, p = p, psi = psi, M = M, Y = Y, Z = Z)
mosy_out <- data.table::CJ(day = 1:tmax, state = c('M', 'Y', 'Z'), value = NaN, species = "mosquito")
mosy_out <- mosy_out[c('M', 'Y', 'Z'), on="state"]
data.table::setkey(mosy_out, day)
human_out <- data.table::CJ(day = 1:tmax, state = c('S', 'I'), value = NaN, species = "human")
human_out <- human_out[c('S', 'I'), on="state"]
data.table::setkey(human_out, day)
while (get_tnow(mod) <= tmax) {
compute_bloodmeal_simple(model = mod)
step_aqua(model = mod)
step_mosquitoes(model = mod)
step_humans(model = mod)
mosy_out[day == get_tnow(mod) & state == 'M', value := mod$mosquito$M]
mosy_out[day == get_tnow(mod) & state == 'Y', value := mod$mosquito$Y]
mosy_out[day == get_tnow(mod) & state == 'Z', value := mod$mosquito$Z]
human_out[day == get_tnow(mod) & state == 'S', value := mod$human$H - mod$human$X]
human_out[day == get_tnow(mod) & state == 'I', value := mod$human$X]
mod$global$tnow <- mod$global$tnow + 1L
}
out <- rbind(mosy_out, human_out)
out[, "run" := as.integer(runid)]
return(out)
})
sto_out <- data.table::rbindlist(sto_out)
ggplot(sto_out) +
geom_line(aes(x = day, y = value, color = state, group = run), alpha = 0.3) +
geom_line(data = det_out, aes(x = day, y = value, color = state)) +
facet_wrap(species ~ state, scales = "free")
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.