library(knitr) opts_chunk$set(cache.path='~/.vignette_cache/') library(knitcitations) cite_options(citation_format = 'pandoc') P <- rprojroot::find_rstudio_root_file
Model after r citep('10.1371/journal.pone.0080091')
.
devtools::load_all() library(ggplot2) library(dplyr) library(doMC) library(tidyr)
Set parameters for the simulator. For now, no real network - just two unconnected patches
#registerDoMC(cores=2) parms = list( beta = 0.004, #contact rate for direct transmission gamma = 0.167, #recovery rate mu = 0, #base mortality rate alpha = 0.111, #disease mortality rate phi = 1.96e-4, #infectiousness of environmental virions eta = 0.14, #degradation rate of environmental virions nu = 0.001, #uptake rate of environmental virion sigma = 0, #virion shedding rate omega = 0.0, #movement rate rho = 0, #contact nonlinearity 0=dens-dependent, 1=freq-dependent lambda = 0, #force of infection from external sources chi = make_net(network_type = "smallworld", network_parms = list(dim = 1, size = 10, nei = 2.33, p = 0.0596, multiple = FALSE, loops = FALSE)), #patch connectivity matrix tau_crit = 3, #critical suveillance time I_crit = 5, #threshold for reporting pi_report = 1,#reporting probability pi_detect = 1 #detection probability ) initial_cond <- matrix(c(100, 0, 0, 0), nrow=10, ncol=4, byrow=TRUE) infected_patches <- c(1,2) initial_cond[, 2] <- 1 initial_cond[, 1] <- initial_cond[, 1] - 1
Run the simulator and plot
#a <- metaflu:::sim_gillespie(as.vector(t(initial_cond)), parms, 1:1000, TRUE) output1 <- mf_sim(init = initial_cond, parameters = parms, times=0:100, n_sims = 25) outputs = output1 %>% group_by(sim, time, patch) %>% summarize(S = sum(population[class=="S"]), pop = sum(population), I = sum(population[class=="I"])) %>% group_by() %>% gather("var", "val", S, pop, I) ggplot(outputs, aes(x=time, y=val, group=paste0(sim,var), col=var)) + facet_wrap(~patch) + geom_line(alpha=0.7) + ylim(0,100)
if(!is.null(parameters[["network_parms"]])) { n_patches = parameters[["network_parms"]][["size"]] } else { n_patches = nrow(parameters[["chi"]]) } #Convert scalar parameters to vectors parms_to_extend <- c("lambda", "tau_crit", "I_crit", "pi_report", "pi_detect") for(parm in parms_to_extend) { parameters[[parm]] = rep_len(parameters[[parm]], n_patches) } initial_cond_vec <- as.vector(t(initial_cond)) a <- sim_gillespie(init=initial_cond_vec, parmlist=parameters, times=0:100, progress=FALSE)
write.bibtex(file="references.bib")
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.