---
title: "Epidemic Risk"
subtitle: "Superspreading on policy, decision making and control measures"
output: rmarkdown::html_vignette
bibliography: references.json
link-citations: true
vignette: >
  %\VignetteIndexEntry{Epidemic Risk}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

::: {.alert .alert-info}
**New to the {superspreading} package?**

See the [Get started](https://epiverse-trace.github.io/superspreading/articles/superspreading.html) vignette for an introduction to the concepts discussed below and some of the functions included in the R package.
:::

During the early stages of an infectious disease outbreak information on transmission
clusters and heterogeneity in individual-level transmission can provide insight into the occurrence of superspreading events and from this
the probability an outbreak will cause an epidemic. This has implications for policy
decisions on how to bring disease transmission 
under control (e.g. $R$ < 1) and increase the probability that an outbreak goes extinct.

A recent example of where individual-level transmission and superspreading were analysed
to understand disease transmission was by the [Scientific Advisory Group for Emergencies 
(SAGE)](https://www.gov.uk/government/organisations/scientific-advisory-group-for-emergencies) 
for the [UK's COVID-19 response](https://www.gov.uk/government/publications/analysis-of-sars-cov-2-transmission-clusters-and-superspreading-events-3-june-2020). 

This vignette explores the applications of the functions included in the {superspreading} package and their visualisation in understanding and informing decision-makers for a variety of outbreak scenarios.

```{r setup}
library(superspreading)
library(ggplot2)
library(ggtext)
```

Early in the COVID-19 epidemic, a number of studies including @kucharskiEarlyDynamicsTransmission2020 investigated
whether data on SARS-CoV-2 incidence had evidence of superspreading events. These events are often estimated as some measure of overdispersion. This was in light of evidence that other corona virus
outbreaks (Severe acute respiratory syndrome (SARS) and Middle East respiratory syndrome (MERS))
had superspreading events. The detection of variability in individual-transmission
is important in understanding or predicting the likelihood of transmission establishing
when imported into a new location (e.g. country). 

First we set up the parameters to use for calculating the probability of an
epidemic. In this case we are varying the heterogeneity of transmission (`k`)
and number of initial infections (`a`).

```{r, prep-dispersion-plot}
epidemic_params <- expand.grid(
  R = 2.35,
  R_lw = 1.5,
  R_up = 3.2,
  k = seq(0, 1, 0.1),
  num_init_infect = seq(0, 10, 1)
)
```

Then we calculate the probability of an epidemic for each parameter combination.
The results are combined with the the parameters.

```{r, calc-prob-endemic}
# results are transposed to pivot to long rather than wide data
prob_epidemic <- t(apply(epidemic_params, 1, function(x) {
  median <- probability_epidemic(
    R = x[["R"]],
    k = x[["k"]],
    num_init_infect = x[["num_init_infect"]]
  )
  lower <- probability_epidemic(
    R = x[["R_lw"]],
    k = x[["k"]],
    num_init_infect = x[["num_init_infect"]]
  )
  upper <- probability_epidemic(
    R = x[["R_up"]],
    k = x[["k"]],
    num_init_infect = x[["num_init_infect"]]
  )
  c(prob_epidemic = median, prob_epidemic_lw = lower, prob_epidemic_up = upper)
}))
epidemic_params <- cbind(epidemic_params, prob_epidemic)
```

To visualise the influence of transmission variation (`k`) we subset the results
to only include those with a single initial infection seeding transmission (`a = 1`).

```{r, subset-prob-epidemic}
# subset data for a single initial infection
homogeneity <- subset(epidemic_params, num_init_infect == 1)
```

```{r, plot-dispersion, fig.cap="The probability that an initial infection (introduction) will cause a sustained outbreak (transmission chain). The dispersion of the individual-level transmission is plotted on the x-axis and probability of outbreak -- calculated using `probability_epidemic()` -- is on the y-axis. This plot is reproduced from @kucharskiEarlyDynamicsTransmission2020 figure 3A.", fig.width = 8, fig.height = 5}
# plot probability of epidemic across dispersion
ggplot(data = homogeneity) +
  geom_ribbon(
    mapping = aes(
      x = k,
      ymin = prob_epidemic_lw,
      ymax = prob_epidemic_up
    ),
    fill = "grey70"
  ) +
  geom_line(mapping = aes(x = k, y = prob_epidemic)) +
  geom_vline(
    mapping = aes(xintercept = 0.2),
    linetype = "dashed"
  ) +
  annotate(geom = "text", x = 0.15, y = 0.75, label = "SARS") +
  geom_vline(
    mapping = aes(xintercept = 0.4),
    linetype = "dashed"
  ) +
  annotate(geom = "text", x = 0.45, y = 0.75, label = "MERS") +
  scale_y_continuous(
    name = "Probability of large outbreak",
    limits = c(0, 1)
  ) +
  scale_x_continuous(name = "Extent of homogeneity in transmission") +
  theme_bw()
```

The degree of variability in transmission increases as $k$ decreases [@lloyd-smithSuperspreadingEffectIndividual2005]. 
So the probability of a large outbreak is smaller for smaller values of $k$, meaning that if COVID-19
is more similar to SARS than MERS it will be less likely to establish and cause an 
outbreak if introduced into a newly susceptible population.

However, this probability of establishment and continued transmission will also
depend on the number of initial introductions to that populations. The more
introductions, the higher the chance one will lead to an epidemic.

```{r, prep-introductions-plot}
introductions <- subset(epidemic_params, k == 0.5)
```

```{r, plot-introductions, fig.cap="The probability that an a number of introduction events will cause a sustained outbreak (transmission chain). The number of disease introductions is plotted on the x-axis and probability of outbreak -- calculated using `probability_epidemic()` -- is on the y-axis. This plot is reproduced from Kucharski et al. (2020) figure 3B.", fig.width = 8, fig.height = 5}
# plot probability of epidemic across introductions
ggplot(data = introductions) +
  geom_pointrange(
    mapping = aes(
      x = num_init_infect,
      y = prob_epidemic,
      ymin = prob_epidemic_lw,
      ymax = prob_epidemic_up
    )
  ) +
  scale_y_continuous(
    name = "Probability of large outbreak",
    limits = c(0, 1)
  ) +
  scale_x_continuous(name = "Number of introductions", n.breaks = 6) +
  theme_bw()
```

Different levels of heterogeneity in transmission will produce different probabilities of epidemics.

```{r, plot-introductions-multi-k, fig.cap="The probability that an a number of introduction events will cause a sustained outbreak (transmission chain). The number of disease introductions is plotted on the x-axis and probability of outbreak -- calculated using `probability_epidemic()` -- is on the y-axis. Different values of dispersion are plotted to show the effect of increased transmission variability on an epidemic establishing", fig.width = 8, fig.height = 5}
# plot probability of epidemic across introductions for multiple k
ggplot(data = epidemic_params) +
  geom_point(
    mapping = aes(
      x = num_init_infect,
      y = prob_epidemic,
      colour = k
    )
  ) +
  scale_y_continuous(
    name = "Probability of large outbreak",
    limits = c(0, 1)
  ) +
  labs(colour = "Dispersion (*k*)") +
  scale_x_continuous(name = "Number of introductions", n.breaks = 6) +
  scale_colour_continuous(type = "viridis") +
  theme_bw() +
  theme(legend.title = element_markdown())
```

In the above plot we drop the uncertainty around each point and assume a known value
of $R$ in order to more clearly show the pattern.

Applications that make it easy to explore this functionality and compare between existing
parameter estimates of offspring distributions are also useful. An example of this is the
[Shiny app developed by the Centre for Mathematical Modelling of Infectious diseases](https://cmmid.github.io/visualisations/new-outbreak-probability) at the London
School of Hygiene and Tropical Medicine, for comparing SARS-like and MERS-like scenarios, 
as well as random mixing.

Conversely to the probability of an epidemics, the probability that an outbreak will 
go extinct (i.e. transmission will subside), can also be plotted for different values of 
$R$.

```{r, prep-exinction-plot}
extinction_params <- expand.grid(
  R = seq(0, 5, 0.1),
  k = c(0.01, 0.1, 0.5, 1, 4, Inf),
  num_init_infect = 1
)
# results are transposed to pivot to long rather than wide data
prob_extinct <- apply(extinction_params, 1, function(x) {
  median <- probability_extinct(
    R = x[["R"]],
    k = x[["k"]],
    num_init_infect = x[["num_init_infect"]]
  )
  median
})
extinction_params <- cbind(extinction_params, prob_extinct)
```

```{r, plot-extinction, fig.cap="The probability that an infectious disease will go extinct for a given value of $R$ and $k$. This is calculated using `probability_extinct()` function. This plot is reproduced from @lloyd-smithSuperspreadingEffectIndividual2005 figure 2B.", fig.width = 8, fig.height = 5}
# plot probability of extinction across R for multiple k
ggplot(data = extinction_params) +
  geom_point(
    mapping = aes(
      x = R,
      y = prob_extinct,
      colour = factor(k)
    )
  ) +
  scale_y_continuous(
    name = "Probability of extinction",
    limits = c(0, 1)
  ) +
  labs(colour = "Dispersion (*k*)") +
  scale_x_continuous(
    name = "Reproductive number (*R*)",
    n.breaks = 6
  ) +
  scale_colour_viridis_d() +
  theme_bw() +
  theme(
    axis.title.x = element_markdown(),
    legend.title = element_markdown()
  )
```

## Superspreading in decision making

One of the first tasks in an outbreak is to establish whether estimates of individual-level
transmission variability have been reported, either for this outbreak or a previous outbreak
of the same pathogen. Alternatively, if it is an outbreak of a novel pathogen, parameters
for similar pathogens can be searched for. 

This is accomplished by another Epiverse-TRACE R package: [{epiparameter}](https://epiverse-trace.github.io/epiparameter/index.html). This package
contains a library of epidemiological parameters included offspring distribution parameter estimates and uncertainty obtained from published literature. This aids the collation and establishes a baseline understanding of transmission. Similar to the collation of transmission clusters for [COVID-19](https://www.gov.uk/government/publications/analysis-of-sars-cov-2-transmission-clusters-and-superspreading-events-3-june-2020).

Equipped with these parameter estimates, metrics that provide an intuitive understanding of transmission patterns can be tabulated.
For example the proportion of transmission events in a given cluster size can be calculated using the {superspreading} function `proportion_cluster_size()`.

This calculates the proportion of new cases that originated within a transmission event of a
given size. In other words, what proportion of all transmission events were part of secondary
case clusters (i.e. from the same primary case) of at least, for example, five cases. This 
calculation is useful to inform backwards contact tracing efforts. 

In this example we look at clusters of at least five, ten and twenty-five secondary cases, at
different numbers of initial infections and for two reproduction numbers to see how this affects
cluster sizes.

```{r}
# For R = 0.8
proportion_cluster_size(
  R = 0.8,
  k = seq(0.1, 1, 0.1),
  cluster_size = c(5, 10, 25)
)

# For R = 3
proportion_cluster_size(
  R = 3,
  k = seq(0.1, 1, 0.1),
  cluster_size = c(5, 10, 25)
)
```

These results show that as the level of heterogeneity in individual-level
transmission increases, a larger percentage of cases come from larger cluster sizes,
and that large clusters can be produced when $R$ is higher even with low levels of 
transmission variation.

It indicates whether preventing gatherings of a certain size can reduce the epidemic
by preventing potential superspreading events.

When data on the settings of transmission is known, priority can be given to restricting 
particular types of gatherings to most effectively bring down the reproduction number.

## Controls on transmission

In their work on identifying the prevalence of overdispersion in individual-level transmission, @lloyd-smithSuperspreadingEffectIndividual2005 also showed the effect of control measures on the probability of containment. They defined containment as the size of a transmission chain not reaching 100 cases. They assumed the implementation of control measures at the individual or population level. {superspreading} provides `probability_contain()` to calculate the probability that an outbreak will go extinct before reaching a threshold size.

```{r, prep-containment-plot}
contain_params <- expand.grid(
  R = 3, k = c(0.1, 0.5, 1, Inf), num_init_infect = 1, control = seq(0, 1, 0.05)
)
prob_contain <- apply(contain_params, 1, function(x) {
  probability_contain(
    R = x[["R"]],
    k = x[["k"]],
    num_init_infect = x[["num_init_infect"]],
    pop_control = x[["control"]]
  )
})
contain_params <- cbind(contain_params, prob_contain)
```

```{r, plot-containment, fig.cap="The probability that an outbreak will be contained (i.e. not exceed 100 cases) for a variety of population-level control measures. The probability of containment is calculated using `probability_contain()`. This plot is reproduced from Lloyd-Smith et al. (2005) figure 3C.", fig.width = 8, fig.height = 5}
# plot probability of epidemic across introductions for multiple k
ggplot(data = contain_params) +
  geom_point(
    mapping = aes(
      x = control,
      y = prob_contain,
      colour = factor(k)
    )
  ) +
  scale_y_continuous(
    name = "Probability of containment",
    limits = c(0, 1)
  ) +
  scale_x_continuous(name = "Control measures (*c*)", n.breaks = 6) +
  labs(colour = "Dispersion (*k*)") +
  scale_colour_viridis_d() +
  theme_bw() +
  theme(
    axis.title.x = element_markdown(),
    legend.title = element_markdown()
  )
```

## References