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

## ----setup--------------------------------------------------------------------
library(superspreading)
library(ggplot2)
library(purrr)
library(dplyr)

## ----load-offspring-dist------------------------------------------------------
library(epiparameter)
offspring_dists <- epiparameter_db(
  epi_name = "offspring distribution"
)

## -----------------------------------------------------------------------------
diseases <- vapply(offspring_dists, `[[`, FUN.VALUE = character(1), "disease")
offspring_dists <- offspring_dists[!duplicated(diseases)]
diseases <- diseases[!duplicated(diseases)]

offspring_dist_params <- vapply(
  offspring_dists, get_parameters, FUN.VALUE = numeric(2)
)

offspring_dist_params <- data.frame(
  disease = diseases,
  t(offspring_dist_params)
)
offspring_dist_params

## -----------------------------------------------------------------------------
offspring_dist_params$t20 <- do.call(
  rbind,
  apply(
    offspring_dist_params,
    MARGIN = 1,
    function(x) {
      proportion_transmission(
        R = as.numeric(x[2]), k = as.numeric(x[3]), percent_transmission = 0.2,
        method = "t_20", format_prop = FALSE
      )
    }
  )
)[, 3]

offspring_dist_params$p80 <- do.call(
  rbind,
  apply(
    offspring_dist_params,
    MARGIN = 1,
    function(x) {
      proportion_transmission(
        R = as.numeric(x[2]), k = as.numeric(x[3]), percent_transmission = 0.8,
        method = "p_80", format_prop = FALSE
      )
    }
  )
)[, 3]
offspring_dist_params

## -----------------------------------------------------------------------------
# nolint start for `:::`
get_infectious_curve <- function(R, k) {
  # upper limit of x when y = 0
  upper_u <- superspreading:::solve_for_u(prop = 0, R = R, k = k)
  upper_u <- round(upper_u)
  u_seq <- superspreading:::lseq(from = 1e-5, to = upper_u, length.out = 500)
  res <- lapply(u_seq, function(upper) {
    integrate(
      function(u) u * superspreading:::fvx(x = u, R, k),
      lower = 0, upper = upper
    )$value / R
  })
  expected_v_more_than_x <- 1 - unlist(res)
  proportion_more_than_x <- 1 - superspreading:::pgammaRk(u_seq, R = R, k = k)
  data.frame(
    exp_t = expected_v_more_than_x,
    prop_i = proportion_more_than_x
  )
}
# nolint end

## ----fig.width=8--------------------------------------------------------------
infect_curve <- map(offspring_dist_params %>%
      group_split(disease), function(x) {
        get_infectious_curve(R = x$mean, k = x$dispersion) %>%
          mutate(
            disease = x$disease,
            R = x$mean, k = x$dispersion
          )
      })

infect_curve <- do.call(rbind, infect_curve)

ggplot(
  data = infect_curve,
  aes(x = prop_i, y = exp_t, colour = disease)
) +
  geom_line() +
  geom_abline(slope = 1) +
  geom_vline(xintercept = 0.2, linetype = 2, alpha = 0.2) +
  scale_x_continuous(breaks = seq(0, 1, by = 0.2)) +
  theme_classic() +
  theme(
    aspect.ratio = 1,
    legend.position = "right"
  ) +
  annotate(
    geom = "text",
    angle = 45,
    size = 2.5,
    x = 0.5,
    y = 0.5,
    vjust = 1.5,
    label = "Homogeneous population"
  ) +
  labs(
    x = "Proportion of infectious cases (ranked)",
    y = "Expected proportion of transmission",
    colour = ""
  ) +
  coord_cartesian(expand = FALSE)

## ----fig.width=8--------------------------------------------------------------
k_seq <- superspreading:::lseq(from = 0.01, to = 100, length.out = 1000) # nolint
y <- map_dbl(
  k_seq,
  function(x) {
    proportion_transmission(
      R = 2, k = x, percent_transmission = 0.2,
      method = "t_20", format_prop = FALSE
    )[, 3]
  }
)
prop_t20 <- data.frame(k_seq, y)

ggplot() +
  geom_line(data = prop_t20, mapping = aes(x = k_seq, y = y)) +
  geom_point(
  data = offspring_dist_params,
    mapping = aes(
      x = dispersion,
      y = t20,
      fill = disease
    ),
    shape = 21,
    size = 3
  ) +
  geom_hline(yintercept = c(0.2, 0.8), linetype = 2, alpha = 0.2) +
  theme_classic() +
  theme(
    aspect.ratio = 1
  ) +
  scale_y_continuous(
    name = paste(
      "Proportion of transmission expected from",
      "the most infectious 20% of cases",
      sep = " \n"
    ),
    limits = c(0, 1),
    breaks = seq(0, 1, by = 0.2)
  ) +
  scale_x_log10(name = "Dispersion parameter", expand = c(0, 0)) +
  labs(fill = "")

## -----------------------------------------------------------------------------
# For k = 0.5
proportion_transmission(
  R = 0.1, k = 0.5, percent_transmission = 0.8, method = "t_20"
)
proportion_transmission(
  R = 1, k = 0.5, percent_transmission = 0.8, method = "t_20"
)
proportion_transmission(
  R = 5, k = 0.5, percent_transmission = 0.8, method = "t_20"
)

# For k = 2
proportion_transmission(
  R = 0.1, k = 2, percent_transmission = 0.8, method = "t_20"
)
proportion_transmission(
  R = 1, k = 2, percent_transmission = 0.8, method = "t_20"
)
proportion_transmission(
  R = 5, k = 2, percent_transmission = 0.8, method = "t_20"
)

## -----------------------------------------------------------------------------
# For k = 0.5
proportion_transmission(
  R = 0.1, k = 0.5, percent_transmission = 0.8, method = "p_80"
)
proportion_transmission(
  R = 1, k = 0.5, percent_transmission = 0.8, method = "p_80"
)
proportion_transmission(
  R = 5, k = 0.5, percent_transmission = 0.8, method = "p_80"
)

# For k = 2
proportion_transmission(
  R = 0.1, k = 2, percent_transmission = 0.8, method = "p_80"
)
proportion_transmission(
  R = 1, k = 2, percent_transmission = 0.8, method = "p_80"
)
proportion_transmission(
  R = 5, k = 2, percent_transmission = 0.8, method = "p_80"
)

## -----------------------------------------------------------------------------
# R = 1, k = 0.5
proportion_transmission(
  R = 1, k = 0.5, percent_transmission = 0.2, method = "p_80"
)
proportion_transmission(
  R = 1, k = 0.5, percent_transmission = 0.8, method = "t_20"
)

# R = 3, k = 2
proportion_transmission(
  R = 3, k = 2, percent_transmission = 0.2, method = "p_80"
)
proportion_transmission(
  R = 3, k = 2, percent_transmission = 0.8, method = "t_20"
)

## -----------------------------------------------------------------------------
1 - proportion_transmission(
  R = 1, k = 0.5, percent_transmission = 0.8, method = "p_80",
  format_prop = FALSE
)[, 3]
proportion_transmission(
  R = 1, k = 0.5, percent_transmission = 0.2, method = "t_20"
)

1 - proportion_transmission(
  R = 1, k = 0.5, percent_transmission = 0.2, method = "t_20",
  format_prop = FALSE
)[, 3]
proportion_transmission(
  R = 1, k = 0.5, percent_transmission = 0.8, method = "p_80"
)

## -----------------------------------------------------------------------------
proportion_transmission(
  R = 1, k = Inf, percent_transmission = 0.8, method = "p_80"
)
proportion_transmission(
  R = 1, k = Inf, percent_transmission = 0.8, method = "t_20"
)