## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 8,
  fig.height = 5
)

## ----setup--------------------------------------------------------------------
library(superspreading)
library(ggplot2)
library(scales)

## ----calc-prob-emerge---------------------------------------------------------
R_wild <- seq(0, 1.2, by = 0.01)
R_mutant <- c(1.2, 1.5, 1000)
mutation_rate <- c(10^-1, 10^-3)

params <- expand.grid(
  R_wild = R_wild,
  R_mutant = R_mutant,
  mutation_rate = mutation_rate
)

prob_emerge <- apply(
  params,
  MARGIN = 1,
  FUN = function(x) {
    probability_emergence(
      R_wild = x[["R_wild"]],
      R_mutant = x[["R_mutant"]],
      mutation_rate = x[["mutation_rate"]],
      num_init_infect = 1
    )
  }
)
res <- cbind(params, prob_emerge = prob_emerge)

## ----plot-prob-emerge---------------------------------------------------------
ggplot(data = res) +
  geom_line(
    mapping = aes(
      x = R_wild,
      y = prob_emerge,
      colour = as.factor(mutation_rate),
      linetype = as.factor(R_mutant)
    )
  ) +
  scale_y_continuous(
    name = "Probability of emergence",
    transform = "log",
    breaks = log_breaks(n = 6),
    limits = c(10^-5, 1)
  ) +
  scale_colour_manual(
    name = "Mutation rate",
    values = c("#228B22", "black")
  ) +
  scale_linetype_manual(
    name = "Mutant pathogen R0",
    values = c("dotted", "dashed", "solid")
  ) +
  scale_x_continuous(
    name = expression(R[0] ~ of ~ introduced ~ pathogen),
    limits = c(0, 1.2),
    n.breaks = 7
  ) +
  theme_bw()

## ----calc-prob-emerge-multi-type----------------------------------------------
R_wild <- seq(0, 1.2, by = 0.01)
R_mutant <- 1.5
mutation_rate <- c(10^-1)
num_mutants <- 0:4

params <- expand.grid(
  R_wild = R_wild,
  R_mutant = R_mutant,
  mutation_rate = mutation_rate,
  num_mutants = num_mutants
)

prob_emerge <- apply(
  params,
  MARGIN = 1,
  FUN = function(x) {
    probability_emergence(
      R_wild = x[["R_wild"]],
      R_mutant = c(rep(x[["R_wild"]], x[["num_mutants"]]), x[["R_mutant"]]),
      mutation_rate = x[["mutation_rate"]],
      num_init_infect = 1
    )
  }
)
res <- cbind(params, prob_emerge = prob_emerge)

## ----plot-prob-emerge-multi-type----------------------------------------------
ggplot(data = res) +
  geom_line(
    mapping = aes(
      x = R_wild,
      y = prob_emerge,
      colour = as.factor(num_mutants)
    )
  ) +
  scale_y_continuous(
    name = "Probability of emergence",
    transform = "log",
    breaks = log_breaks(n = 6),
    limits = c(10^-5, 1)
  ) +
  scale_colour_brewer(
    name = "Number of \nintermediate mutants",
    palette = "Set1"
  ) +
  scale_x_continuous(
    name = expression(R[0] ~ of ~ introduced ~ pathogen),
    limits = c(0, 1.2),
    n.breaks = 7
  ) +
  theme_bw()

## ----calc-prob-emerge-multi-init-infect---------------------------------------
R_wild <- seq(0, 1.2, by = 0.01)
R_mutant <- 1.5
mutation_rate <- 10^-1
num_mutants <- 1
num_init_infect <- seq(1, 9, by = 2)

params <- expand.grid(
  R_wild = R_wild,
  R_mutant = R_mutant,
  mutation_rate = mutation_rate,
  num_mutants = num_mutants,
  num_init_infect = num_init_infect
)

prob_emerge <- apply(
  params,
  MARGIN = 1,
  FUN = function(x) {
    probability_emergence(
      R_wild = x[["R_wild"]],
      R_mutant = c(rep(x[["R_wild"]], x[["num_mutants"]]), x[["R_mutant"]]),
      mutation_rate = x[["mutation_rate"]],
      num_init_infect = x[["num_init_infect"]]
    )
  }
)
res <- cbind(params, prob_emerge = prob_emerge)

## ----plot-prob-emerge-multi-init-infect---------------------------------------
ggplot(data = res) +
  geom_line(
    mapping = aes(
      x = R_wild,
      y = prob_emerge,
      colour = as.factor(num_init_infect)
    )
  ) +
  scale_y_continuous(
    name = "Probability of emergence",
    transform = "log",
    breaks = log_breaks(n = 6),
    limits = c(10^-5, 1)
  ) +
  scale_colour_brewer(
    name = "Number of introductions",
    palette = "Set1"
  ) +
  scale_x_continuous(
    name = expression(R[0] ~ of ~ introduced ~ pathogen),
    limits = c(0, 1.2),
    n.breaks = 7
  ) +
  theme_bw()