## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

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

## ----calc-R-zika--------------------------------------------------------------
infect_duration <- exp(seq(log(0.01), log(10), length.out = 100))
prob_transmission <- exp(seq(log(0.01), log(1), length.out = 100))
params <- expand.grid(
  infect_duration = infect_duration,
  prob_transmission = prob_transmission
)
R <- t(apply(
  params,
  MARGIN = 1,
  function(x) {
    calc_network_R(
      mean_num_contact = 14.1,
      sd_num_contact = 69.6,
      infect_duration = x[["infect_duration"]],
      prob_transmission = x[["prob_transmission"]],
      age_range = c(16, 74)
    )
  }
))
res <- cbind(params, R)
res <- reshape(
  data = res,
  varying = c("R", "R_net"),
  v.names = "R",
  direction = "long",
  times = c("R", "R_net"),
  timevar = "group"
)

## ----plot-zika-r, fig.cap="The reproduction number using the unadjusted and adjusted calculation -- calculated using `calc_network_R()` -- with mean duration of infection on the x-axis and transmission probability per sexual partner on the y-axis. The line shows the points that $R_0$ is equal to one. Both axes are plotted on a natural log scale. This plot is similar to Figure 1 from @yakobLowRiskSexuallytransmitted2016, but is plotted as a heatmap and without annotation.", fig.width = 8, fig.height = 5----
ggplot(data = res) +
  geom_tile(
    mapping = aes(
      x = infect_duration,
      y = prob_transmission,
      fill = R
    )
  ) +
  geom_contour(
    mapping = aes(
      x = infect_duration,
      y = prob_transmission,
      z = R
    ),
    breaks = 1,  # Set the contour line at R = 1
    colour = "black"
  ) +
  scale_x_continuous(
    name = "Mean duration of infectiousness",
    trans = "log", breaks = breaks_log()
  ) +
  scale_y_continuous(
    name = "Zika virus transmission probability per sexual partners",
    trans = "log", breaks = breaks_log()
  ) +
  facet_wrap(
    vars(group),
    labeller = as_labeller(c(R = "R", R_net = "Adjusted R"))
  ) +
  scale_fill_viridis_c(
    trans = "log",
    breaks = breaks_log(n = 8)(min(res$R):max(res$R)),
    labels = label_comma()
  ) +
  theme_bw()

## ----calc-r-mpox--------------------------------------------------------------
beta <- seq(0.001, 1, length.out = 1000)
duration_years <- 21 / 365
res <- lapply(
  beta,
  calc_network_R,
  mean_num_contact = 10,
  sd_num_contact = 50,
  infect_duration = duration_years,
  age_range = c(18, 44)
)
res <- do.call(rbind, res)
res <- as.data.frame(cbind(beta, res))
res <- reshape(
  data = res,
  varying = c("R", "R_net"),
  v.names = "R",
  direction = "long",
  times = c("R", "R_net"),
  timevar = "group"
)

## ----plot-mpox, fig.cap="The reproduction number using the unadjusted and adjusted calculation -- calculated using `calc_network_R()` -- with secondary attack rate on the x-axis and reproduction number ($R_0$) on the y-axis. This plot is similar to Figure 2A from @endoHeavytailedSexualContact2022.", fig.width = 8, fig.height = 5----
ggplot(data = res) +
  geom_line(mapping = aes(x = beta, y = R, colour = group)) +
  geom_hline(mapping = aes(yintercept = 1)) +
  scale_y_continuous(
    name = "Reproduction Number (R)",
    trans = "log",
    breaks = breaks_log(),
    labels = label_comma()
  ) +
  scale_x_continuous(name = "Secondary Attack Rate (SAR)") +
  scale_colour_brewer(
    name = "Reproduction Number",
    labels = c("R", "Adjusted R"),
    palette = "Set1"
  ) +
  theme_bw() +
  theme(legend.position = "top")