# Modelling uncertainty in R₀

Epidemic final size calculations are sensitive to input data such as the \(R_0\) of the infection. Such values can often be uncertain in the early stages of an outbreak. This uncertainty can be included in final size calculations by running `final_size()` for values drawn from a distribution, and summarising the outcomes.

New to finalsize? It may help to read earlier vignettes first!

## Use case

The infection parameter (\(R_0\)) is uncertain. We want to know how much variation this could cause in the estimated final size of the epidemic.

### What we have

1. In addition to the \(R_0\), demography data, social contact data, and data on the distribution of susceptibility among age groups;
2. A measure of the error in the estimated \(R_0\) of the infection.

### What we assume

1. An SIR epidemic, and the complete partitioning of individuals into different demographic and infection risk groups.

Code

    library(finalsize)

    # load necessary packages
    library(data.table)
    library(ggplot2)

## Getting \(R_0\), contact and demography data, and susceptibility

This example uses social contact data from the POLYMOD project to estimate the final size of an epidemic in the U.K. These data are provided with the `socialmixr` package.

These data are handled just as in the Basic usage vignette. This example also considers an infection with an \(R_0\) of 2.0.

Code

    # get UK polymod data
    polymod <- socialmixr::polymod
    contact_data <- socialmixr::contact_matrix(
      polymod,
      countries = "United Kingdom",
      age.limits = c(0, 5, 18, 40, 65),
      symmetric = TRUE
    )

    # get the contact matrix and demography data
    contact_matrix <- t(contact_data$matrix)
    demography_vector <- contact_data$demography$population

    # scale the contact matrix so the largest eigenvalue is 1.0
    contact_matrix <- contact_matrix / max(eigen(contact_matrix)$values)

    # divide each row of the contact matrix by the corresponding demography
    contact_matrix <- contact_matrix / demography_vector

    n_demo_grps <- length(demography_vector)

Code

    # mean R0 is 2.0
    r0_mean <- 2.0

For simplicity, this example considers a scenario in which susceptibility to infection does not vary.

Code

    # susceptibility is uniform
    susc_uniform <- matrix(
      data = 1,
      nrow = n_demo_grps,
      ncol = 1L
    )

    # p_susceptibility is uniform
    p_susc_uniform <- susc_uniform

## Running `final_size` over \(R_0\) samples

The basic reproduction number \(R_0\) of an infection might be uncertain in the early stages of an epidemic. This uncertainty can be modelled by running `final_size()` multiple times for the same contact, demography, and susceptibility data, while sampling \(R_0\) values from a distribution.

This example assumes that the \(R_0\) estimate, and the uncertainty around that estimate, is provided as the mean and standard deviation of a normal distribution.

This example considers a normal distribution \(N(\mu = 2.0, \sigma = 0.1)\), for an \(R_0\) of 2.0. We can draw 1,000 \(R_0\) samples from this distribution and run `final_size()` on the contact data and demography data for each sample.

This is quick, as `finalsize` is an Rcpp package with a C++ backend.

Code

    # create an R0 samples vector
    r0_samples <- rnorm(n = 1000, mean = r0_mean, sd = 0.1)

Code

    # run final size on each sample with the same data
    final_size_data <- Map(
      r0_samples, seq_along(r0_samples),
      f = function(r0, i) {
        # the i-th final size estimate
        fs <- final_size(
          r0 = r0,
          contact_matrix = contact_matrix,
          demography_vector = demography_vector,
          susceptibility = susc_uniform,
          p_susceptibility = p_susc_uniform
        )

        fs$replicate <- i
        fs$r0_estimate <- r0
        fs
      }
    )

    # combine data
    final_size_data <- rbindlist(final_size_data)

    # order age groups
    final_size_data$demo_grp <- factor(
      final_size_data$demo_grp,
      levels = contact_data$demography$age.group
    )

### Visualise uncertainty in final size

Code

    ggplot(final_size_data) +
      stat_summary(
        aes(
          demo_grp, p_infected
        ),
        fun = mean,
        fun.min = function(x) {
          quantile(x, 0.05)
        },
        fun.max = function(x) {
          quantile(x, 0.95)
        }
      ) +
      scale_y_continuous(
        labels = scales::percent,
        limits = c(0.25, 1)
      ) +
      theme_classic() +
      theme(
        legend.position = "top",
        legend.key.height = unit(2, "mm"),
        legend.title = ggtext::element_markdown(
          vjust = 1
        )
      ) +
      coord_cartesian(
        expand = TRUE
      ) +
      labs(
        x = "Age group",
        y = "% Infected"
      )

Figure 1: Estimated ranges of the final size of a hypothetical SIR epidemic in age groups of the U.K. population, when the \(R_0\) is estimated to be 2.0, with a standard deviation around this estimate of 0.1. In this example, relatively low uncertainty in \(R_0\) estimates can also lead to uncertainty in the estimated final size of the epidemic. Points represent means, while ranges extend between the 5th and 95th percentiles.