In ecohealthalliance/metaflu: Functions for Avian Influenza Metapopulation Simulation

library(knitr)
opts_chunk$set(cache.path='~/.vignette_cache/', echo = FALSE, message=FALSE, warning = FALSE, autodep=TRUE, cache = TRUE)
library(knitcitations)
cite_options(citation_format = 'pandoc')
P <- rprojroot::find_rstudio_root_file
source(P("R/summarizing_functions.R"))

Model after r citep('10.1371/journal.pone.0080091').

library(metaflu)
library(ggplot2)
library(dplyr)
library(doMC)
library(tidyr)

get_outputs <- function(funcname, matrow){
  print(paste("Starting",matrow[1],"by",matrow[2],"simulation."))
  results <- do.call(funcname, list(matrow[1], matrow[2]))
  duration <- get_duration(results)
  lower_d <- quantile(unlist(duration$days_greater), probs = c(0.025))
  upper_d <- quantile(unlist(duration$days_greater), probs = c(0.975))
  median_d <- quantile(unlist(duration$days_greater), probs = c(0.5))
  failure <- proportion_failed(get_failure(results))
  abundance <- get_tot_infections(results)
  lower_a <- quantile(unlist(abundance$total_i), probs = c(0.025))
  upper_a <- quantile(unlist(abundance$total_i), probs = c(0.975))
  median_a <- quantile(unlist(abundance$total_i), probs = c(0.5))
  df <- data.frame(failure, lower_d, median_d, upper_d, lower_a, median_a, upper_a)
  df
}

set.seed(17)
fig2 <- function(farm_size, farm_number){
  initial_cond <- matrix(c(farm_size, 0, 0, 0), nrow=farm_number, ncol=4, byrow=TRUE)
  infected_patches <- sample(seq_len(nrow(initial_cond)), 2)
  initial_cond[infected_patches, 2] <- 1
  initial_cond[infected_patches, 1] <- initial_cond[infected_patches, 1] - 1

  fig2parms = list(
    beta = 0.004,   #contact rate for direct transmission
    gamma = 0.167,  #recovery rate
    mu = 0,         #base mortality rate
    alpha = 0.1111,      #disease mortality rate
    phi = 0,  #infectiousness of environmental virions
    eta = 0,     #degradation rate of environmental virions
    nu =  0.00,    #uptake rate of environmental virion
    sigma = 0,      #virion shedding rate
    omega = 0.03,   #movement rate
    rho = 0,        #contact  nonlinearity 0=dens-dependent, 1=freq-dependent
    lambda = 0,     #force of infection from external sources
    tau_crit = 0,   #critical suveillance time
    I_crit = 0,     #threshold for reporting
    pi_report = .0, #reporting probability
    pi_detect = .0, #detection probability
    chi = make_net(network_type = "smallworld",
                   network_parms = list(dim = 1, size = farm_number, nei = 2.33, p = 0.0596, multiple = FALSE, loops = FALSE)),
    stochastic_network = TRUE
  )

  x <- mf_sim(init = initial_cond, parameters = fig2parms, times=0:1000, n_sims = 100)
  x
}

registerDoMC(cores=34)

recreate_data <- function(funcname){
  a <- c(40,50,64,80,100,125,160,200,256,320,400,500,640,800)
  parMat <- cbind(a, rev(a))
  results <- apply(parMat, 1, function(x) get_outputs(funcname,x))
  results
}

fig2results <- bind_rows(recreate_data("fig2"))

saveRDS(fig2results, "fig2results.rds")

#full_recreate_data <- function(funcname){
 # a <- c(5,8,10,16,20,25,32,40,50,64,80,100,125,160,200,256,320,400,500,640,800,1000,1280,1600,2000,3200,4000,6400)
#  parMat <- cbind(a, rev(a))
#  results <- apply(parMat, 1, function(x) get_outputs(funcname, x))
#  results
#}

#fullfig2results <- bind_rows(full_recreate_data("fig2"))

ggplot(data = fig2results) +
  geom_line(aes(x = seq_along(failure), y = failure)) +
  scale_x_continuous(breaks = seq_along(fig2results$failure), name = "Chickens:Farm", labels = c("40:800", "50:640", "64:500", "80:400","100:320", "125:256", "160:200","200:160","256:125","320:100","400:80", "500:64", "640:50", "800:40")) +
  scale_y_continuous(name = "Proportion of Epidemic Failures") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1))

The above figure is a recreation of the top graph of Figure 2, showing the 'probability' that an epidemic will fail using their model; in this context "epidemic failure" is a run where the epidemic does not propagate an originally seeded patch to any other patch. Our model shows roughly the same behavior as the published paper, with epidemic failure becoming less likely as we increase the size of farms.

ggplot(data = fig2results) +
  geom_line(aes(x = seq_along(median_d), y = median_d)) +
  geom_line(aes(x = seq_along(lower_d), y = lower_d), linetype = "dashed") +
  geom_line(aes(x = seq_along(upper_d), y = upper_d), linetype = "dashed") + 
  scale_x_continuous(breaks = seq_along(fig2results$failure), name = "Chickens:Farm", labels = c("40:800", "50:640", "64:500", "80:400","100:320", "125:256", "160:200","200:160","256:125","320:100","400:80", "500:64", "640:50", "800:40")) +
  scale_y_continuous(name = "Median Duration (days)") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1))

This graph is a recreation of the middle graph in Figure 2, and shows the median length of the epidemic in days, along with the 95% bootstrapped interval in dashed lines. This graph is considerably different from the published results, with epidemics that last far shorter and simulations with a much lower 97.5% quantile. Although Figure 2 includes environmental transmission, Figure S8 does not, and still reports a median duration in the 300/500 range for these Chickens:Farm combinations. '

ggplot(data = fig2results) +
  geom_line(aes(x = seq_along(median_a), y = median_a)) +
  geom_line(aes(x = seq_along(lower_a), y = lower_a), linetype = "dashed") +
  geom_line(aes(x = seq_along(upper_a), y = upper_a), linetype = "dashed") + 
  scale_x_continuous(breaks = seq_along(fig2results$failure), name = "Chickens:Farm", labels = c("40:800", "50:640", "64:500", "80:400","100:320", "125:256", "160:200","200:160","256:125","320:100","400:80", "500:64", "640:50", "800:40")) +
  scale_y_continuous(name = "Number of Infections") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1))

The above graph shows the median total number of infections as we alter the number of chickens and number of farms. The fact that it looks a bit like the graph of a binary variable reflects the fact that the results of each run at each combination of chickens/farms either propagate to nearly the entire population (~32000) or barely propagate at all (~5-15 infections).

write.bibtex(file="references.bib")

References

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