In ecohealthalliance/metaflu: Functions for Avian Influenza Metapopulation Simulation

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")

References

ecohealthalliance/metaflu documentation built on May 15, 2019, 7:56 p.m.