Outbreak analytics with R: All in One View

Content from Outbreak analytics pipelines

Last updated on 2024-09-24 | Edit this page

Estimated time: 12 minutes

Overview

Questions

Why use R packages for Outbreak analytics?
What can we do to analyse our outbreak data?
How can I start doing Outbreak Analytics with R?

Objectives

Explain our vision on the need for outbreak analytics R packages.
Share our strategy to create R packages into an outbreak analytics pipeline.
Define our plan to start your learning path in outbreak analytics with R.

Prerequisites

This episode requires you to be familiar with:

Data science : Basic programming with R.

Epidemic theory : Reproduction number.

Why to use R packages for Outbreak analytics?

Outbreaks appear with different diseases and in different contexts, but what all of them have in common are the key public health questions (Cori et al. 2017).

Is the epidemic going to take off? Is it under control? How much effort will be needed to control it? We can answer them by quantifying the transmissibility of the disease. The most used parameter for this is the reproduction number ($R$), the average number of secondary infections caused by a typical primary case in the population of interest (Prism, 2016). We can intuitively interpret it as:

if $R>1$, the epidemic is likely to grow,
if $R<1$, the epidemic is likely to decline.

We can estimate the reproduction number by initially using two data inputs: the incidence of reported cases and the generation time distribution. But to calculate it, we must apply the appropriate mathematical models written in code with the required computational methods. That is not enough! Following good practices, the code we write should be peer-reviewed and contain internal tests to double-check that we are getting the estimates we expect. Imagine rewriting all of it during a health emergency!

In R, the fundamental unit of shareable code is the package. A package bundles together code, data, documentation, and tests and is easy to share with others (Wickham and Bryan, 2023). We, as epidemiologists, can contribute to their collaborative maintenance as a community to perform less error-prone data analysis pipelines.

Questions to think about

Remember your last experience with outbreak data and reflect on these questions:

What data sources did you need to understand the outbreak?
How did you get access to that data?
Is that analysis pipeline you followed reusable for the next response?

Reflect on your experiences.

Example: Quantify transmission

The EpiNow2 package provides a three-step solution to quantify the transmissibility. Let’s see how to do this with a minimal example. First, load the package:

R

library(EpiNow2)

First, get your case data

Case incidence data must be stored in a data frame with the observed number of cases per day. We can read an example from the package:

R

example_confirmed

OUTPUT

           date confirm
         <Date>   <num>
  1: 2020-02-22      14
  2: 2020-02-23      62
  3: 2020-02-24      53
  4: 2020-02-25      97
  5: 2020-02-26      93
 ---
126: 2020-06-26     296
127: 2020-06-27     255
128: 2020-06-28     175
129: 2020-06-29     174
130: 2020-06-30     126

Then, set the generation time

Not all primary cases have the same probability of generating a secondary case. The onset and cessation of infectiousness may occur gradually. For EpiNow2, we can specify it as a probability distribution with mean, standard deviation sd, and maximum value max:

R

generation_time <- dist_spec(
  mean = 3.6,
  sd = 3.1,
  max = 20,
  distribution = "lognormal"
)

Let’s calculate the reproduction number!

In the epinow() function we can add:

the reported_cases data frame,
the generation_time delay distribution, and
the computation stan parameters for this calculation:

R

epinow_estimates <- epinow(
  # cases
  reported_cases = example_confirmed[1:60],
  # delays
  generation_time = generation_time_opts(generation_time),
  # computation
  stan = stan_opts(
    cores = 4, samples = 1000, chains = 3,
    control = list(adapt_delta = 0.99)
  )
)

OUTPUT

WARN [2024-09-24 01:16:09] epinow: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
Running the chains for more iterations may help. See
https://mc-stan.org/misc/warnings.html#bulk-ess -
WARN [2024-09-24 01:16:10] epinow: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
Running the chains for more iterations may help. See
https://mc-stan.org/misc/warnings.html#tail-ess -

As an output, we get the time-varying (or effective) reproduction number, as well as the cases by date of report and date of infection:

R

base::plot(epinow_estimates)

Is this $Rt$ estimation biased?

Review Gostic et al., 2020 about what additional adjustments this estimation requires to avoid false precision in $Rt$.

The problem!

However, quantifying the transmissibility during a real-life outbreak response is more challenging than this example!

Usually, we receive outbreak data in non-standard formats, requiring specific steps and taking the most time to prepare usable data inputs. Some of them are:

Read delay distributions from the literature
Read and clean case data
Validate your line list
Describe case data

And this is not the end. After quantifying transmissibility we need to answer more key public health questions like: What is the attack rate we expect? What would be the impact of a given intervention? We can use the reproduction number and other outputs as new inputs for complementary tasks. For example:

Estimate severity
Create short-term forecast
Simulate transmission scenarios
Compare interventions

So, all these tasks can be interconnected in a pipeline:

What can we do?

Our strategy is gradually incorporating specialised R packages into our traditional analysis pipeline. These packages should fill the gaps in these epidemiology-specific tasks in response to outbreaks.

Epiverse-TRACE’s aim is to provide a software ecosystem for outbreak analytics. We support the development of software pieces, make the existing ones interoperable for the user experience, and stimulate a community of practice.

How can I start?

Our plan for these tutorials is to introduce key solutions from packages in all the tasks before and after the Quantify transmission task, plus the required theory concepts to interpret modelling outputs and make rigorous conclusions.

In the first set of episodes, you will learn how to optimise the reading of delay distributions and cleaning of case data to input them into the Quantify transmission task. These preliminary tasks are the Early tasks. These include packages like {readepi}, cleanepi, linelist, {epiparameter}, and {episoap}.
Then, we will get deeper into the packages and required theory to Quantify transmission and perform more real-time analysis tasks next to it. These are the Middle tasks. This includes EpiNow2, cfr, {epichains}, and {superspreading}.
Lastly, we will use Quantify transmission data outputs to compare it to other indicators and simulate epidemic scenarios as part of the Late tasks. This includes finalsize, {epidemics}, and {scenarios}.

Key Points

Our vision is to have pipelines of R packages for outbreak analytics.
Our strategy is to create interconnected tasks to get relevant outputs for public health questions.
We plan to introduce package solutions and theory bits for each of the tasks in the outbreak analytics pipeline.

Content from Read delays

Last updated on 2024-09-24 | Edit this page

Estimated time: 12 minutes

Overview

Questions

How to get delay distributions from a systematic review?
How to connect reused delays with my existing analysis pipeline?
When should delays be reused from a systematic review?

Objectives

Get delays from a systematic review with {epiparameter}.
Get statistical summaries and distribution parameters of delay distributions.
Use distribution functions from delay distributions.
Convert a continuous to a discrete delay distribution.

Prerequisites

This episode requires you to be familiar with:

Data science : Basic programming with R.

Epidemic theory : Epidemiological parameters. Time periods.

Introduction

The natural history of an infectious disease shows that its development has a regularity from stage to stage. The time periods from an infectious disease inform about the timing of transmission and interventions.

Definition of key time periods. From Xiang et al, 2021

Definitions

Look at the glossary for the definitions of all the time periods of the figure above!

However, early in an epidemic, modelling efforts can be delayed by the lack of a centralized resource that summarises input parameters for the disease of interest (Nash et al., 2023). Projects like {epiparameter} and {epireview} are building online catalogues following systematic review protocols that can help build models faster for coming outbreaks and epidemics from known pathogens and unknown ones related to known families of viruses.

To exemplify how to use {epiparameter} in your analysis pipeline, our goal in this episode will be to replace the generation_time input that we can use for EpiNow2::epinow().

R

epinow_estimates <- epinow(
  # cases
  reported_cases = example_confirmed[1:60],
  # delays
  generation_time = generation_time_opts(generation_time),
  # computation
  stan = stan_opts(
    cores = 4, samples = 1000, chains = 3,
    control = list(adapt_delta = 0.99)
  )
)

To do this replacement, instead of plug-in numeric values to EpiNow2::dist_spec() to manually specify the delay distribution parameters, we are going to collect them from the library of epidemiological parameters provided by {epiparameter}:

R

generation_time <- dist_spec(
  mean = 3.6,
  sd = 3.1,
  max = 20,
  distribution = "lognormal"
)

Let’s explore how we can access this and other time delays using {epiparameter}. We’ll use the pipe %>% to connect some of their functions, so let’s also call to the tidyverse package:

R

library(epiparameter)
library(EpiNow2)
library(tidyverse)

Find a Generation time

The generation time, jointly with the $R$, can inform about the speed of spread and its feasibility of control. Given a $R>1$, with a shorter generation time, cases can appear more quickly.

Video from the MRC Centre for Global Infectious Disease Analysis, Ep 76. Science In Context - Epi Parameter Review Group with Dr Anne Cori (27-07-2023) at https://youtu.be/VvpYHhFDIjI?si=XiUyjmSV1gKNdrrL

In calculating the effective reproduction number ($R_{t}$), the generation time distribution is often approximated by the serial interval distribution. This frequent approximation is because it is easier to observe and measure the onset of symptoms than the onset of infectiousness.

A schematic of the relationship of different time periods of transmission between an infector and an infectee in a transmission pair. Exposure window is defined as the time interval having viral exposure, and transmission window is defined as the time interval for onward transmission with respect to the infection time (Chung Lau et al. 2021).

However, using the serial interval as an approximation of the generation time is primarily valid for diseases in which infectiousness starts after symptom onset (Chung Lau et al. 2021). In cases where infectiousness starts before symptom onset, the serial intervals can have negative values, which is the case of a pre-symptomatic transmission (Nishiura et al. (2020)).

Additionally, even if the generation time and serial interval have the same mean, their variance usually differs, propagating bias to the $R_{t}$ estimation. $R_{t}$ estimates are sensitive not only to the mean generation time but also to the variance and form of the generation interval distribution (Gostic et al., 2020).

From time periods to probability distributions.

When we calculate the serial interval, we see that not all case pairs have the same time length. We will observe this variability for any case pair and individual time period, including the incubation period and infectious period.

Serial intervals of possible case pairs in (a) COVID-19 and (b) MERS-CoV. Pairs represent a presumed infector and their presumed infectee plotted by date of symptom onset (Althobaity et al., 2022).

To summarize these data from individual and pair time periods, we can find the statistical distributions that best fit the data (McFarland et al., 2023).

Fitted serial interval distribution for (a) COVID-19 and (b) MERS-CoV based on reported transmission pairs in Saudi Arabia. We fitted three commonly used distributions, Lognormal, Gamma, and Weibull distributions, respectively (Althobaity et al., 2022).

Statistical distributions are summarized in terms of their summary statistics like the location (mean and percentiles) and spread (variance or standard deviation) of the distribution, or with their distribution parameters that inform about the form (shape and rate/scale) of the distribution. These estimated values can be reported with their uncertainty (95% confidence intervals).

Gamma	mean	shape	rate/scale
MERS-CoV	14.13(13.9–14.7)	6.31(4.88–8.52)	0.43(0.33–0.60)
COVID-19	5.1(5.0–5.5)	2.77(2.09–3.88)	0.53(0.38–0.76)

Weibull	mean	shape	rate/scale
MERS-CoV	14.2(13.3–15.2)	3.07(2.64–3.63)	16.1(15.0–17.1)
COVID-19	5.2(4.6–5.9)	1.74(1.46–2.11)	5.83(5.08–6.67)

Serial interval estimates using Gamma, Weibull, and Log normal distributions. 95% confidence intervals for the shape and scale (logmean and sd for Log normal) parameters are shown in brackets (Althobaity et al., 2022).
Log normal	mean	mean-log	sd-log
MERS-CoV	14.08(13.1–15.2)	2.58(2.50–2.68)	0.44(0.39–0.5)
COVID-19	5.2(4.2–6.5)	1.45(1.31–1.61)	0.63(0.54–0.74)

Serial interval

Assume that COVID-19 and SARS have similar reproduction number values and that the serial interval approximates the generation time.

Given the Serial interval of both infections in the figure below:

Which one would be harder to control?
Why do you conclude that?

Serial interval of novel coronavirus (COVID-19) infections overlaid with a published distribution of SARS. (Nishiura et al, 2020)

Give me a hint

The peak of each curve can inform you about the location of the mean of each distribution. The larger the mean, the larger the serial interval.

Show me the solution

Which one would be harder to control?

COVID-19

Why do you conclude that?

COVID-19 has the lowest mean serial interval. The approximate mean value for the serial interval of COVID-19 is around four days, and SARS is about seven days. Thus, COVID-19 will likely have newer generations in less time than SARS, assuming similar reproduction numbers.

Instructor Note

The objective of the assessment above is to assess the interpretation of a larger or shorter generation time.

Extract epidemiological parameters

First, let’s assume that the data set example_confirmed has COVID-19 observed cases. So, we need to find a reported generation time for COVID-19 or any other useful parameter for this aim.

Let’s start by looking at how many parameters we have in the epidemiological distributions database (epidist_db) for the disease named covid-19:

R

epiparameter::epidist_db(
  disease = "covid"
)

OUTPUT

Returning 27 results that match the criteria (22 are parameterised).
Use subset to filter by entry variables or single_epidist to return a single entry.
To retrieve the short citation for each use the 'get_citation' function

OUTPUT

List of <epidist> objects
  Number of entries in library: 27
  Number of studies in library: 10
  Number of diseases: 1
  Number of delay distributions: 27
  Number of offspring distributions: 0

From the {epiparameter} package, we can use the epidist_db() function to ask for any disease and also for a specific epidemiological distribution (epi_dist).

Let’s ask now how many parameters we have in the epidemiological distributions database (epidist_db) with the generation time using the string generation:

R

epiparameter::epidist_db(
  epi_dist = "generation"
)

OUTPUT

Returning 1 results that match the criteria (1 are parameterised).
Use subset to filter by entry variables or single_epidist to return a single entry.
To retrieve the short citation for each use the 'get_citation' function

OUTPUT

Disease: Influenza
Pathogen: Influenza-A-H1N1
Epi Distribution: generation time
Study: Lessler J, Reich N, Cummings D, New York City Department of Health and
Mental Hygiene Swine Influenza Investigation Team (2009). "Outbreak of
2009 Pandemic Influenza A (H1N1) at a New York City School." _The New
England Journal of Medicine_. doi:10.1056/NEJMoa0906089
<https://doi.org/10.1056/NEJMoa0906089>.
Distribution: weibull
Parameters:
  shape: 2.360
  scale: 3.180

Currently, in the library of epidemiological parameters, we have one generation time entry for Influenza. Considering the abovementioned considerations, we can look at the serial intervals for COVID-19.

R

epiparameter::epidist_db(
  disease = "COVID",
  epi_dist = "serial"
)

OUTPUT

Returning 4 results that match the criteria (3 are parameterised).
Use subset to filter by entry variables or single_epidist to return a single entry.
To retrieve the short citation for each use the 'get_citation' function

OUTPUT

[[1]]
Disease: COVID-19
Pathogen: SARS-CoV-2
Epi Distribution: serial interval
Study: Alene M, Yismaw L, Assemie M, Ketema D, Gietaneh W, Birhan T (2021).
"Serial interval and incubation period of COVID-19: a systematic review
and meta-analysis." _BMC Infectious Diseases_.
doi:10.1186/s12879-021-05950-x
<https://doi.org/10.1186/s12879-021-05950-x>.
Parameters: <no parameters>

[[2]]
Disease: COVID-19
Pathogen: SARS-CoV-2
Epi Distribution: serial interval
Study: Nishiura H, Linton N, Akhmetzhanov A (2020). "Serial interval of novel
coronavirus (COVID-19) infections." _International Journal of
Infectious Diseases_. doi:10.1016/j.ijid.2020.02.060
<https://doi.org/10.1016/j.ijid.2020.02.060>.
Distribution: lnorm
Parameters:
  meanlog: 1.386
  sdlog: 0.568

[[3]]
Disease: COVID-19
Pathogen: SARS-CoV-2
Epi Distribution: serial interval
Study: Nishiura H, Linton N, Akhmetzhanov A (2020). "Serial interval of novel
coronavirus (COVID-19) infections." _International Journal of
Infectious Diseases_. doi:10.1016/j.ijid.2020.02.060
<https://doi.org/10.1016/j.ijid.2020.02.060>.
Distribution: weibull
Parameters:
  shape: 2.203
  scale: 5.420

[[4]]
Disease: COVID-19
Pathogen: SARS-CoV-2
Epi Distribution: serial interval
Study: Yang L, Dai J, Zhao J, Wang Y, Deng P, Wang J (2020). "Estimation of
incubation period and serial interval of COVID-19: analysis of 178
cases and 131 transmission chains in Hubei province, China."
_Epidemiology and Infection_. doi:10.1017/S0950268820001338
<https://doi.org/10.1017/S0950268820001338>.
Distribution: norm
Parameters:
  mean: 4.600
  sd: 4.400

attr(,"class")
[1] "multi_epidist"

CASE-INSENSITIVE

epidist_db is case-insensitive. This means that you can use strings with letters in upper or lower case indistinctly.

We get more than one epidemiological delay. To summarize this view and get the column names from the underlying parameter dataset, we can add the epiparameter::list_distributions() function to the previous code using the pipe %>%:

R

epiparameter::epidist_db(
  disease = "covid",
  epi_dist = "serial"
) %>%
  epiparameter::list_distributions()

OUTPUT

Returning 4 results that match the criteria (3 are parameterised).
Use subset to filter by entry variables or single_epidist to return a single entry.
To retrieve the short citation for each use the 'get_citation' function

OUTPUT

   disease epi_distribution prob_distribution       author year
1 COVID-19  serial interval              <NA> Muluneh .... 2021
2 COVID-19  serial interval             lnorm Hiroshi .... 2020
3 COVID-19  serial interval           weibull Hiroshi .... 2020
4 COVID-19  serial interval              norm Lin Yang.... 2020

Ebola’s incubation periods

Take 5 minutes:

How many delay distributions are for the Ebola disease?
How many delay distributions are for the incubation period of Ebola?
Explore the library and find the disease with the delay distribution of your interest! Do you recognize the paper?

Give me a hint

The {epiparameter} combo of epidist_db() plus list_distributions() list all the entries by:

disease,
epidemiological distribution,
the type of the probability distribution,
author of the study, and
year of study.

Show me the solution

R

# 16 delays distributions
epiparameter::epidist_db(
  disease = "ebola"
)

# 5 delay distributions are for the incubation period
epiparameter::epidist_db(
  disease = "ebola",
  epi_dist = "incubation"
)

Now, from the output of epiparameter::epidist_db(), What is an offspring distribution?

Select a single distribution

The epiparameter::epidist_db() function works as a filtering or subset function. Let’s use the author argument to filter Hiroshi Nishiura parameters:

R

epiparameter::epidist_db(
  disease = "covid",
  epi_dist = "serial",
  author = "Hiroshi"
) %>%
  epiparameter::list_distributions()

OUTPUT

Returning 2 results that match the criteria (2 are parameterised).
Use subset to filter by entry variables or single_epidist to return a single entry.
To retrieve the short citation for each use the 'get_citation' function

OUTPUT

   disease epi_distribution prob_distribution       author year
1 COVID-19  serial interval             lnorm Hiroshi .... 2020
2 COVID-19  serial interval           weibull Hiroshi .... 2020

We still get more than one epidemiological parameter. We can set the single_epidist argument to TRUE to only one:

R

epiparameter::epidist_db(
  disease = "covid",
  epi_dist = "serial",
  author = "Hiroshi",
  single_epidist = TRUE
)

OUTPUT

Using Nishiura H, Linton N, Akhmetzhanov A (2020). "Serial interval of novel
coronavirus (COVID-19) infections." _International Journal of
Infectious Diseases_. doi:10.1016/j.ijid.2020.02.060
<https://doi.org/10.1016/j.ijid.2020.02.060>..
To retrieve the short citation use the 'get_citation' function

OUTPUT

Disease: COVID-19
Pathogen: SARS-CoV-2
Epi Distribution: serial interval
Study: Nishiura H, Linton N, Akhmetzhanov A (2020). "Serial interval of novel
coronavirus (COVID-19) infections." _International Journal of
Infectious Diseases_. doi:10.1016/j.ijid.2020.02.060
<https://doi.org/10.1016/j.ijid.2020.02.060>.
Distribution: lnorm
Parameters:
  meanlog: 1.386
  sdlog: 0.568

How does `single_epidist` works?

Looking at the help documentation for ?epiparameter::epidist_db():

If multiple entries match the arguments supplied and single_epidist = TRUE,
Then, the parameterised ⁠<epidist>⁠ with the largest sample size will be returned.
If multiple entries are equal after this sorting, the first entry will be returned.

What does a parametrised <epidist> is? Look at ?is_parameterised.

Now, we have an epidemiological parameter we can reuse! We can replace the numbers we plug into EpiNow2::dist_spec().

Let’s assign this <epidist> class object to the covid_serialint object.

R

covid_serialint <-
  epiparameter::epidist_db(
    disease = "covid",
    epi_dist = "serial",
    author = "Nishiura",
    single_epidist = TRUE
  )

OUTPUT

Using Nishiura H, Linton N, Akhmetzhanov A (2020). "Serial interval of novel
coronavirus (COVID-19) infections." _International Journal of
Infectious Diseases_. doi:10.1016/j.ijid.2020.02.060
<https://doi.org/10.1016/j.ijid.2020.02.060>..
To retrieve the short citation use the 'get_citation' function

R

covid_serialint

OUTPUT

Disease: COVID-19
Pathogen: SARS-CoV-2
Epi Distribution: serial interval
Study: Nishiura H, Linton N, Akhmetzhanov A (2020). "Serial interval of novel
coronavirus (COVID-19) infections." _International Journal of
Infectious Diseases_. doi:10.1016/j.ijid.2020.02.060
<https://doi.org/10.1016/j.ijid.2020.02.060>.
Distribution: lnorm
Parameters:
  meanlog: 1.386
  sdlog: 0.568

Ebola’s incubation period

Take 2 minutes:

What type of distribution has the incubation period of Ebola with the highest sample size?
How would you access to the sample size of the other studies in the <multi_epidist> class object?

Give me a hint

The {epiparameter} combo of epidist_db() plus list_distributions() list all the entries by:

disease,
epidemiological distribution,
the type of the probability distribution,
author of the study, and
year of study.

This is a <multi_epidist> class object:

R

epiparameter::epidist_db(
  disease = "ebola",
  epi_dist = "incubation"
)

Show me the solution

R

# the distribution with the highest sample size has a gamma distribution
epiparameter::epidist_db(
  disease = "ebola",
  epi_dist = "incubation",
  single_epidist = TRUE
)

To access the sample_size, review an issue reported in the GitHub repository of the {epiparameter} package.

Extract the summary statistics

We can get the mean and standard deviation (sd) from this <epidist> diving into the summary_stats object:

R

# get the mean
covid_serialint$summary_stats$mean

OUTPUT

[1] 4.7

How to get the `sd` and other nested elements?

Take 1 minute to:

Get the sd of the epidemiological distribution.
Find the sample_size used in the study.
Explore all the other nested elements within the <epidist> object.

Share about:

What elements do you find useful for your analysis?
What other elements would you like to see in this object? How?

Give me a hint

Use the $ operator plus the tab keyboard button to explore them as an expandable list:

R

covid_serialint$

Use the str() to display the structure of the <epidist> R object.

Show me the solution

R

# get the sd
covid_serialint$summary_stats$sd

# get the sample_size
covid_serialint$metadata$sample_size

An interesting element is the method_assess nested entry, which refers to the methods used by the study authors to assess for bias while estimating the serial interval distribution.

R

covid_serialint$method_assess

OUTPUT

$censored
[1] TRUE

$right_truncated
[1] TRUE

$phase_bias_adjusted
[1] FALSE

We will explore these concepts at the end!

Ebola’s severity parameter

A severity parameter like the duration of hospitalization could add to the information needed about the bed capacity in response to an outbreak (Cori et al., 2017).

For Ebola:

what is a reported point estimate and uncertainty of the mean duration of health-care and case isolation?

Give me a hint

An informative delay measures the time from symptom onset to recovery or death.

Show me the solution

R

# one way to get the list of all the available parameters
epidist_db(disease = "all") %>%
  list_distributions() %>%
  as_tibble() %>%
  distinct(epi_distribution)

ebola_severity <- epidist_db(
  disease = "ebola",
  epi_dist = "onset to discharge"
)

# point estimate
ebola_severity$summary_stats$mean
# 95% confidence intervals
ebola_severity$summary_stats$mean_ci
# limits of the confidence intervals
ebola_severity$summary_stats$mean_ci_limits

Continuous distributions

The following output has four entries with different content in the probability distribution (prob_distribution) column:

R

distribution <-
  epiparameter::epidist_db(
    disease = "covid",
    epi_dist = "serial"
  )

OUTPUT

Returning 4 results that match the criteria (3 are parameterised).
Use subset to filter by entry variables or single_epidist to return a single entry.
To retrieve the short citation for each use the 'get_citation' function

R

distribution %>%
  list_distributions()

OUTPUT

   disease epi_distribution prob_distribution       author year
1 COVID-19  serial interval              <NA> Muluneh .... 2021
2 COVID-19  serial interval             lnorm Hiroshi .... 2020
3 COVID-19  serial interval           weibull Hiroshi .... 2020
4 COVID-19  serial interval              norm Lin Yang.... 2020

Entries with a missing value (<NA>) in the prob_distribution column are non-parameterised entries. They have summary statistics but no probability distribution. Compare these two outputs:

R

distribution[[1]]$summary_stats
distribution[[1]]$prob_dist

Parameterised entries have an Inference method

As detailed in ?is_parameterised, a parameterised distribution is the entry that has a probability distribution associated with it provided by an inference_method as shown in metadata:

R

distribution[[1]]$metadata$inference_method
distribution[[2]]$metadata$inference_method
distribution[[4]]$metadata$inference_method

In the epiparameter::list_distributions() output, we can also find different types of probability distributions (e.g., Log-normal, Weibull, Normal).

R

distribution %>%
  list_distributions()

OUTPUT

   disease epi_distribution prob_distribution       author year
1 COVID-19  serial interval              <NA> Muluneh .... 2021
2 COVID-19  serial interval             lnorm Hiroshi .... 2020
3 COVID-19  serial interval           weibull Hiroshi .... 2020
4 COVID-19  serial interval              norm Lin Yang.... 2020

In {epiparameter}, you will mostly find continuous distributions like these. You can visualize any of them with the plot() function and access to:

the Probability Density Function (PDF) and
the Cumulative Distribution Function (CDF).

R

plot(distribution[[2]])

With the day_range argument, you can change the length or number of days in the x axis. Explore what it look like:

R

plot(distribution[[2]], day_range = 0:20)

The distribution Zoo

Explore this shinyapp called The Distribution Zoo!

Follow these steps to reproduce the form of the covid_serialint distribution:

Access to https://ben18785.shinyapps.io/distribution-zoo/ shinyapp website,
Go to the left panel,
Keep the Category of distribution: Continuous Univariate,
Select a new Type of distribution: Log-Normal,
Move the sliders, i.e. the graphical control element that allows you to adjust a value by moving a handle along a horizontal track or bar to the covid_serialint parameters.

Replicate these with the distribution object and all its list elements: 2, 3, and 4. Explore how the shape of a distribution changes when its parameters change.

Share about:

What other features of the website do you find helpful?

Distribution functions

In R, all the statistical distributions have functions to access the:

Probability Density function (PDF),
Cumulative Distribution function (CDF),
Quantile function, and
Random values from the given distribution.

Functions for the Normal distribution

If you need it, read in detail about the R probability functions for the normal distribution, each of its definitions and identify in which part of a distribution they are located!

If you look at ?stats::Distributions, each type of distribution has a unique set of functions. However, {epiparameter} gives you the same four functions to access each of the values above for any <epidist> object you want!

R

# plot this to have a visual reference
plot(covid_serialint, day_range = 0:20)

R

# the density value at quantile value of 10 (days)
density(covid_serialint, at = 10)

OUTPUT

[1] 0.01911607

R

# the cumulative probability at quantile value of 10 (days)
cdf(covid_serialint, q = 10)

OUTPUT

[1] 0.9466605

R

# the quantile value (day) at a cumulative probability of 60%
quantile(covid_serialint, p = 0.6)

OUTPUT

[1] 4.618906

R

# generate 10 random values (days) given
# the distribution family and its parameters
generate(covid_serialint, times = 10)

OUTPUT

 [1] 2.271797 5.917676 6.223395 6.321504 1.785465 4.854480 2.087282 3.985876
 [9] 3.777108 4.779382

Instructor Note

Access to the reference documentation (Help files) for these functions is accessible with the three double-colon notation: epiparameter:::

?epiparameter:::density.epidist()
?epiparameter:::cdf.epidist()
?epiparameter:::quantile.epidist()
?epiparameter:::generate.epidist()

Window for contact tracing and the Serial interval

The serial interval is important in the optimization of contact tracing since it provides a time window for the containment of a disease spread (Fine, 2003). Depending on the serial interval, we can evaluate the need to expand the number of days pre-onset to consider in the contact tracing to include more backwards contacts (Davis et al., 2020).

With the COVID-19 serial interval (covid_serialint) calculate:

How much more of the backward cases could be captured if the contact tracing method considered contacts up to 6 days pre-onset compared to 2 days pre-onset?

Give me a hint

In Figure 5 from the R probability functions for the normal distribution, the shadowed section represents a cumulative probability of 0.997 for the quantile value at x = 2.

Show me the solution

R

plot(covid_serialint)

R

cdf(covid_serialint, q = 2)
cdf(covid_serialint, q = 6)

Given the COVID-19 serial interval:

A contact tracing method considering contacts up to 2 days pre-onset will capture around 11.1% of backward cases.
If this period is extended to 6 days pre-onset, this could include 76.2% of backward contacts.

What if

If we exchange the question between days and cumulative probability to:

When considering secondary cases, how many days following the symptom onset of primary cases can we expect 55% of symptom onset to occur?

R

quantile(covid_serialint, p = 0.55)

An interpretation could be:

The 55% percent of the symptom onset of secondary cases will happen after 4.2 days after the symptom onset of primary cases.

Discretize a continuous distribution

We are getting closer to the end! EpiNow2::dist_spec() still needs a maximum value (max).

One way to do this is to get the quantile value for the distribution’s 99.9th percentile or 0.999 cumulative probability. For this, we need access to the set of distribution functions for our <epidist> object.

We can use the set of distribution functions for a continuous distribution (as above). However, these values will be continuous numbers. We can discretize the continuous distribution stored in our <epidist> object to get discrete values from a continuous distribution.

When we epiparameter::discretise() the continuous distribution we get a discrete(-ized) distribution:

R

covid_serialint_discrete <-
  epiparameter::discretise(covid_serialint)

covid_serialint_discrete

OUTPUT

Disease: COVID-19
Pathogen: SARS-CoV-2
Epi Distribution: serial interval
Study: Nishiura H, Linton N, Akhmetzhanov A (2020). "Serial interval of novel
coronavirus (COVID-19) infections." _International Journal of
Infectious Diseases_. doi:10.1016/j.ijid.2020.02.060
<https://doi.org/10.1016/j.ijid.2020.02.060>.
Distribution: discrete lnorm
Parameters:
  meanlog: 1.386
  sdlog: 0.568

We identify this change in the Distribution: output line of the <epidist> object. Take a double check to this line:

Distribution: discrete lnorm

While for a continuous distribution, we plot the Probability Density Function (PDF), for a discrete distribution, we plot the Probability Mass Function (PMF):

R

# continuous
plot(covid_serialint)

# discrete
plot(covid_serialint_discrete)

To finally get a max value, let’s access the quantile value of the 99.9th percentile or 0.999 probability of the distribution with the prob_dist$q notation, similarly to how we access the summary_stats values.

R

covid_serialint_discrete_max <-
  covid_serialint_discrete$prob_dist$q(p = 0.999)

Lenght of quarantine and Incubation period

The incubation period distribution is a useful delay to assess the length of active monitoring or quarantine (Lauer et al., 2020). Similarly, delays from symptom onset to recovery (or death) will determine the required duration of health-care and case isolation (Cori et al., 2017).

Calculate:

Within what exact time frame do 99% of individuals who develop COVID-19 symptoms exhibit them after infection?

Give me a hint

What delay distribution measures the time between infection and the onset of symptoms?

The probability function for <epidist> discrete distributions differ from the continuous ones!

R

# plot to have a visual reference
plot(covid_serialint_discrete, day_range = 0:20)

# density value at quantile value 10 (day)
covid_serialint_discrete$prob_dist$d(10)

# cumulative probability at quantile value 10 (day)
covid_serialint_discrete$prob_dist$cdf(10)

# In what quantile value (days) do we have the 60% cumulative probability?
covid_serialint_discrete$prob_dist$q(0.6)

# generate random values
covid_serialint_discrete$prob_dist$r(10)

Show me the solution

R

covid_incubation <-
  epiparameter::epidist_db(
    disease = "covid",
    epi_dist = "incubation",
    single_epidist = TRUE
  )

covid_incubation_discrete <- epiparameter::discretise(covid_incubation)

covid_incubation_discrete$prob_dist$q(0.99)

99% of those who develop COVID-19 symptoms will do so within 16 days of infection.

Now, Is this result expected in epidemiological terms?

how to create a distribution plot?

From a maximum value with $prob_dist$q(), we can create a sequence of quantile values as a numeric vector and map density values for each:

R

# create a discrete distribution visualization
# from a maximum value from the distribution
covid_serialint_discrete$prob_dist$q(0.999) %>%
  # generate quantile values
  # as a sequence for each natural number
  seq(1L, to = ., by = 1L) %>%
  # coerce numeric vector to data frame
  as_tibble_col(column_name = "quantile_values") %>%
  mutate(
    # map density values
    # for each quantile in the density function
    density_values =
      covid_serialint_discrete$prob_dist$d(quantile_values)
  ) %>%
  # create plot
  ggplot(
    aes(
      x = quantile_values,
      y = density_values
    )
  ) +
  geom_col()

Plug-in `{epiparameter}` to `{EpiNow2}`

Now we can plug everything into the EpiNow2::dist_spec() function!

R

serial_interval_covid <-
  dist_spec(
    mean = covid_serialint$summary_stats$mean,
    sd = covid_serialint$summary_stats$sd,
    max = covid_serialint_discrete_max,
    distribution = "lognormal"
  )

serial_interval_covid

OUTPUT


  Fixed distribution with PMF [0.18 0.11 0.08 0.066 0.057 0.05 0.045 0.041 0.037 0.034 0.032 0.03 0.028 0.027 0.025 0.024 0.023 0.022 0.021 0.02 0.019 0.019 0.018]

Warning

Using the serial interval instead of the generation time is an alternative that can propagate bias in your estimates, even more so in diseases with reported pre-symptomatic transmission.

Let’s replace the generation_time input we used for EpiNow2::epinow().

R

epinow_estimates <- epinow(
  # cases
  reported_cases = example_confirmed[1:60],
  # delays
  generation_time = generation_time_opts(serial_interval_covid),
  # computation
  stan = stan_opts(
    cores = 4, samples = 1000, chains = 3,
    control = list(adapt_delta = 0.99)
  )
)

base::plot(epinow_estimates)

Ebola’s effective reproduction number

Download and read the Ebola dataset:

Reuse one epidemiological parameter to estimate the effective reproduction number for the Ebola dataset.
Why did you choose that parameter?

Give me a hint

To calculate the $R_t$, we need:

data set with confirmed cases per day and
one key delay distribution

Key functions we applied in this episode are:

epidist_db()
list_distributions()
discretise()
probability functions for continuous and discrete distributions

Show me the solution

R

# read data
# e.g.: if path to file is data/raw-data/ebola_cases.csv then:
ebola_confirmed <-
  read_csv(here::here("data", "raw-data", "ebola_cases.csv"))

# list distributions
epidist_db(disease = "ebola") %>%
  list_distributions()

# subset one distribution
ebola_serial <- epidist_db(
  disease = "ebola",
  epi_dist = "serial",
  single_epidist = TRUE
)

ebola_serial_discrete <- discretise(ebola_serial)

serial_interval_ebola <-
  dist_spec(
    mean = ebola_serial$summary_stats$mean,
    sd = ebola_serial$summary_stats$sd,
    max = ebola_serial_discrete$prob_dist$q(p = 0.999),
    distribution = "gamma"
  )

# name of the type of distribution
# only for the discretised distribution
ebola_serial_discrete$prob_dist$name

epinow_estimates <- epinow(
  # cases
  reported_cases = ebola_confirmed,
  # delays
  generation_time = generation_time_opts(serial_interval_ebola),
  # computation
  stan = stan_opts(
    cores = 4, samples = 1000, chains = 3,
    control = list(adapt_delta = 0.99)
  )
)

plot(epinow_estimates)

Adjusting for reporting delays

Estimating $R_t$ requires data on the daily number of new infections. Due to lags in the development of detectable viral loads, symptom onset, seeking care, and reporting, these numbers are not readily available. All observations reflect transmission events from some time in the past. In other words, if $d$ is the delay from infection to observation, then observations at time $t$ inform $R_{t−d}$, not $R_t$. (Gostic et al., 2020)

Timeline for chain of disease reporting, the Netherlands. Lab, laboratory; PHA, public health authority. From Marinović et al., 2015

The delay distribution could be inferred jointly with the underlying times of infection or estimated as the sum of the incubation period distribution and the distribution of delays from symptom onset to observation from line list data (reporting delay).

For EpiNow2, we can specify these two complementary delay distributions in the delays argument.

Rt is a measure of transmission at time t. Observations after time t must be adjusted. ICU, intensive care unit. From Gostic et al., 2020

Reuse an Incubation period for COVID-19

Use {epiparameter} to:

Find an incubation period for COVID-19.
Add our last epinow() code chunk using the delays argument and the delay_opts() helper function.

Give me a hint

The delays argument and the delay_opts() helper function are analogous to the generation_time argument and the generation_time_opts() helper function.

R

epinow_estimates <- epinow(
  # cases
  reported_cases = example_confirmed[1:60],
  # delays
  generation_time = generation_time_opts(serial_interval_covid),
  delays = delay_opts(incubation_time_covid),
  # computation
  stan = stan_opts(
    cores = 4, samples = 1000, chains = 3,
    control = list(adapt_delta = 0.99)
  )
)

Show me the solution

R

covid_incubation <- epiparameter::epidist_db(
  disease = "covid",
  epi_dist = "incubation",
  author = "Natalie",
  single_epidist = TRUE
)

covid_incubation

covid_incubation_discrete <- epiparameter::discretise(covid_incubation)

incubation_time_covid <- dist_spec(
  mean = covid_incubation$summary_stats$mean,
  sd = covid_incubation$summary_stats$sd,
  max = covid_incubation_discrete$prob_dist$q(p = 0.999),
  distribution = "lognormal"
)

epinow_estimates <- epinow(
  # cases
  reported_cases = example_confirmed[1:60],
  # delays
  generation_time = generation_time_opts(serial_interval_covid),
  delays = delay_opts(incubation_time_covid),
  # computation
  stan = stan_opts(
    cores = 4, samples = 1000, chains = 3,
    control = list(adapt_delta = 0.99)
  )
)

base::plot(epinow_estimates)

How much has it changed?

After adding the incubation period, discuss:

Does the retrospective trend of forecast change?
Has the uncertainty changed?
How would you explain or interpret any of these changes?

Ebola’s effective reproduction number was adjusted by reporting delays

Using the same Ebola dataset:

Reuse one additional epidemiological parameter for the delays argument in EpiNow2::epinow().
Estimate the effective reproduction number using EpiNow2::epinow().
Why did you choose that parameter?

Give me a hint

We can use two complementary delay distributions to estimate the $R_t$ at time $t$.

Show me the solution

R

# read data
# e.g.: if path to file is data/raw-data/ebola_cases.csv then:
ebola_confirmed <-
  read_csv(here::here("data", "raw-data", "ebola_cases.csv"))

# list distributions
epidist_db(disease = "ebola") %>%
  list_distributions()

# subset one distribution for the generation time
ebola_serial <- epidist_db(
  disease = "ebola",
  epi_dist = "serial",
  single_epidist = TRUE
)

ebola_serial_discrete <- discretise(ebola_serial)

serial_interval_ebola <-
  dist_spec(
    mean = ebola_serial$summary_stats$mean,
    sd = ebola_serial$summary_stats$sd,
    max = ebola_serial_discrete$prob_dist$q(p = 0.999),
    distribution = "gamma"
  )

# subset one distribution for delay of the incubation period
ebola_incubation <- epidist_db(
  disease = "ebola",
  epi_dist = "incubation",
  single_epidist = TRUE
)

ebola_incubation_discrete <- discretise(ebola_incubation)

incubation_period_ebola <-
  dist_spec(
    mean = ebola_incubation$summary_stats$mean,
    sd = ebola_incubation$summary_stats$sd,
    max = ebola_incubation_discrete$prob_dist$q(p = 0.999),
    distribution = "gamma"
  )

epinow_estimates <- epinow(
  # cases
  reported_cases = ebola_confirmed,
  # delays
  generation_time = generation_time_opts(serial_interval_ebola),
  delays = delay_opts(incubation_period_ebola),
  # computation
  stan = stan_opts(
    cores = 4, samples = 1000, chains = 3,
    control = list(adapt_delta = 0.99)
  )
)

plot(epinow_estimates)

Extract parameters

Use the influenza_england_1978_school from the outbreaks package to calculate the effective reproduction number.

Give me a hint

How to get the mean and standard deviation from a generation time with median and quantiles as summary statistics?

Look at how to extract parameters from {epiparameter} vignette on parameter extraction and conversion

Show me the solution

R

# What parameters are available for Influenza?
epidist_db(disease = "influenza") %>%
  list_distributions() %>%
  as_tibble() %>%
  count(epi_distribution)

influenza_generation <-
  epidist_db(
    disease = "influenza",
    epi_dist = "generation"
  )

influenza_generation_discrete <-
  discretise(influenza_generation)

# problem
# the summary statistics do not have mean and sd
influenza_generation$summary_stats
influenza_generation$summary_stats$median
influenza_generation$summary_stats$quantiles

# solution
# extract parameters from percentiles
influenza_extracted <- extract_param(
  type = "percentiles",
  values = c(influenza_generation$summary_stats$quantiles[1],
             influenza_generation$summary_stats$quantiles[2]),
  distribution = "lnorm",
  percentiles = c(0.05, 0.95)
)

influenza_extracted

generation_time_influenza <-
  dist_spec(
    mean = influenza_extracted[1],
    sd = influenza_extracted[2],
    max = influenza_generation_discrete$prob_dist$q(p = 0.999),
    distribution = "lognormal"
  )

influenza_cleaned <-
  outbreaks::influenza_england_1978_school %>%
  select(date, confirm = in_bed)

epinow_estimates <- epinow(
  # cases
  reported_cases = influenza_cleaned,
  # delays
  generation_time = generation_time_opts(generation_time_influenza),
  # computation
  stan = stan_opts(
    cores = 4, samples = 1000, chains = 3,
    control = list(adapt_delta = 0.99)
  )
)

plot(epinow_estimates)

When to reuse? When to estimate?

In the early stage of an outbreak, we can rely on reusing parameters for known pathogens to unknown ones, like for the Disease X, a pathogen currently unknown to cause human disease and potentially cause a serious international epidemic (WHO, 2018).

But when data from lines list paired with contact tracing is available, we can estimate the key delay distributions that best fit our data. These will help us to inform, compare and update any previous estimate about questions like:

How long should contacts be followed?
What is the required duration of contact tracing?
How long should cases be isolated to reduce transmission?

However, the methods to accurately estimate delays like the generation interval from contact tracing data involve adjusting for biases like censoring, right truncation and epidemic phase bias. (Gostic et al., 2020)

We can identify what entries in the {epiparameter} library assessed for these biases in their methodology with the method_assess nested entry:

R

covid_serialint$method_assess

OUTPUT

$censored
[1] TRUE

$right_truncated
[1] TRUE

$phase_bias_adjusted
[1] FALSE

How to estimate delay distributions for Disease X?

Refer to this excellent tutorial on estimating the serial interval and incubation period of Disease X accounting for censoring using Bayesian inference with packages like rstan and coarseDataTools.

Tutorial in English: https://rpubs.com/tracelac/diseaseX
Tutorial en Español: https://epiverse-trace.github.io/epimodelac/EnfermedadX.html

How to assess for pre-symptomatic infection?

The lengths of the Serial interval and Incubation period determine the type of disease transmission.

The relationship between the incubation period and serial interval. From Nishiura 2020

Estimating the proportion of pre-symptomatic infections, or the extent to which infectiousness precedes symptom onset will determine the effectiveness of contact tracing and the feasibility of controlling an outbreak (Fraser et al., 2004 and Hellewell et al., 2020).

Parameter estimates. Plausible ranges for the key parameters R0 and θ (read the main text for sources) for four viral infections of public concern are shown as shaded regions. The size of the shaded area reflects the uncertainties in the parameter estimates. Fraser et al., 2004

Meta-analysis on the proportion of pre-symptomatic and asymptomatic transmission in SARS-CoV-2 found limitations of the evidence given high heterogeneity and high risk of selection and information bias between studies (Buitrago-Garcia et al., 2022). This is a call to action to improve the Outbreak Analytic pipelines to use and reuse in the early phase of an outbreak.

What type of transmission?

Compare the serial interval and incubation period of Influenza and MERS:

What type of transmission has Influenza?
What type of transmission has MERS?
Do these results correlate with the available evidence?

Give me a hint

For types of transmission, we refer to infections with symptomatic or pre-symptomatic transmission.

Key functions:

epidist_db()
epidist$summary_stats$

Show me the solution

In this solution we use purrr::pluck() to extract elements within the summary_stats object which is of class list.

R

# pre-symptomatic transmission
epidist_db(
  disease = "influenza",
  epi_dist = "incubation",
  single_epidist = TRUE
) %>%
  pluck("summary_stats") %>%
  pluck("mean")

epidist_db(
  disease = "influenza",
  epi_dist = "serial",
  single_epidist = TRUE
) %>%
  pluck("summary_stats") %>%
  pluck("mean")

# symptomatic transmission
epidist_db(
  disease = "mers",
  epi_dist = "incubation",
  single_epidist = TRUE
) %>%
  pluck("summary_stats") %>%
  pluck("median")

epidist_db(
  disease = "mers",
  epi_dist = "serial",
  single_epidist = TRUE
) %>%
  pluck("summary_stats") %>%
  pluck("mean")

R

# pre-symptomatic transmission
epidist_db(
  disease = "covid",
  epi_dist = "incubation",
  author = "Stephen",
  single_epidist = TRUE
) %>%
  pluck("summary_stats") %>%
  pluck("mean")

epidist_db(
  disease = "covid",
  epi_dist = "serial",
  author = "Nishiura",
  single_epidist = TRUE
) %>%
  pluck("summary_stats") %>%
  pluck("mean")

# symptomatic transmission
epidist_db(
  disease = "ebola",
  epi_dist = "incubation",
  single_epidist = TRUE
) %>%
  pluck("summary_stats") %>%
  pluck("mean")

epidist_db(
  disease = "ebola",
  epi_dist = "serial",
  single_epidist = TRUE
) %>%
  pluck("summary_stats") %>%
  pluck("mean")

Key Points

Use {epiparameter} to access the systematic review catalogue of epidemiological delay distributions.
Use epidist_db() to select single delay distributions.
Use list_distributions() for an overview of multiple delay distributions.
Use discretise() to convert continuous to discrete delay distributions.
Use {epiparameter} probability functions for any delay distributions.

Content from Quantifying transmission

Last updated on 2024-09-24 | Edit this page

Estimated time: 30 minutes

Overview

Questions

How can I estimate key transmission metrics from a time series of case data?
How can I quantify geographical heterogeneity in these metrics?

Objectives

Learn how to estimate transmission metrics from a time series of case data using the R package EpiNow2

Prerequisites

Learners should familiarise themselves with following concept dependencies before working through this tutorial:

Statistics : probability distributions, principle of Bayesian analysis.

Epidemic theory : Effective reproduction number.

Reminder: the Effective Reproduction Number, $R_t$

The basic reproduction number, $R_0$, is the average number of cases caused by one infectious individual in a entirely susceptible population.

But in an ongoing outbreak, the population does not remain entirely susceptible as those that recover from infection are typically immune. Moreover, there can be changes in behaviour or other factors that affect transmission. When we are interested in monitoring changes in transmission we are therefore more interested in the value of the effective reproduction number, $R_t$, the average number of cases caused by one infectious individual in the population at time $t$.

Introduction

Quantifying transmission metrics at the start of an outbreak can give important information on the strength of transmission (reproduction number) and the speed of transmission (growth rate, doubling/halving time). To estimate these key metrics using case data we must account for delays between the date of infections and date of reported cases. In an outbreak situation, data are usually available on reported dates only, therefore we must use estimation methods to account for these delays when trying to understand changes in transmission over time.

In the next tutorials we will focus on how to implement the functions in EpiNow2 to estimate transmission metrics of case data. We will not cover the theoretical background of the models or inference framework, for details on these concepts see the vignette. For more details on the distinction between speed and strength of transmission and implications for control, see Dushoff & Park, 2021.

Bayesian inference

The R package EpiNow2 uses a Bayesian inference framework to estimate reproduction numbers and infection times based on reporting dates.

In Bayesian inference, we use prior knowledge (prior distributions) with data (in a likelihood function) to find the posterior probability.

Posterior probability $\propto$ likelihood $\times$ prior probability

The first step is to load the EpiNow2 package :

R

library(EpiNow2)

Instructor Note

This tutorial illustrates the usage of epinow() to estimate the time-varying reproduction number and infection times. Learners should understand the necessary inputs to the model and the limitations of the model output.

Delay distributions and case data

Case data

To illustrate the functions of EpiNow2 we will use outbreak data of the start of the COVID-19 pandemic from the United Kingdom. The data are available in the R package incidence2.

R

head(incidence2::covidregionaldataUK)

OUTPUT

        date          region region_code cases_new cases_total deaths_new
1 2020-01-30   East Midlands   E12000004        NA          NA         NA
2 2020-01-30 East of England   E12000006        NA          NA         NA
3 2020-01-30         England   E92000001         2           2         NA
4 2020-01-30          London   E12000007        NA          NA         NA
5 2020-01-30      North East   E12000001        NA          NA         NA
6 2020-01-30      North West   E12000002        NA          NA         NA
  deaths_total recovered_new recovered_total hosp_new hosp_total tested_new
1           NA            NA              NA       NA         NA         NA
2           NA            NA              NA       NA         NA         NA
3           NA            NA              NA       NA         NA         NA
4           NA            NA              NA       NA         NA         NA
5           NA            NA              NA       NA         NA         NA
6           NA            NA              NA       NA         NA         NA
  tested_total
1           NA
2           NA
3           NA
4           NA
5           NA
6           NA

To use the data, we must format the data to have two columns:

date : the date (as a date object see ?is.Date()),
confirm : number of confirmed cases on that date.

R

cases <- aggregate(
  cases_new ~ date,
  data = incidence2::covidregionaldataUK[, c("date", "cases_new")],
  FUN = sum
)
colnames(cases) <- c("date", "confirm")

There are case data available for 489 days, but in an outbreak situation it is likely we would only have access to the beginning of this data set. Therefore we assume we only have the first 90 days of this data.

Delay distributions

We assume there are delays from the time of infection until the time a case is reported. We specify these delays as distributions to account for the uncertainty in individual level differences. The delay can consist of multiple types of delays/processes. A typical delay from time of infection to case reporting may consist of :

time from infection to symptom onset (the incubation period) + time from symptom onset to case notification (the reporting time) .

The delay distribution for each of these processes can either estimated from data or obtained from the literature. We can express uncertainty about what the correct parameters of the distributions by assuming the distributions have fixed parameters or whether they have variable parameters. To understand the difference between fixed and variable distributions, let’s consider the incubation period.

Delays and data

The number of delays and type of delay is a flexible input that depends on the data. The examples below highlight how the delays can be specified for different data sources:

Data source	Delay(s)
Time of case report	Incubation period + time from symptom onset to case notification
Time of hospitalisation	Incubation period + time from symptom onset to hospitalisation
Time of symptom onset	Incubation period

Incubation period distribution

The distribution of incubation period can usually be obtained from the literature. The package {epiparameter} contains a library of epidemiological parameters for different diseases obtained from the literature.

We specify a (fixed) gamma distribution with mean $\mu = 4$ and standard deviation $\sigma^2= 2$ (shape = $4$, scale = $1$) using the function dist_spec() as follows:

R

incubation_period_fixed <- dist_spec(
  mean = 4, sd = 2,
  max = 20, distribution = "gamma"
)
incubation_period_fixed

OUTPUT


  Fixed distribution with PMF [0.019 0.12 0.21 0.21 0.17 0.11 0.069 0.039 0.021 0.011 0.0054 0.0026 0.0012 0.00058 0.00026 0.00012 5.3e-05 2.3e-05 1e-05 4.3e-06]

The argument max is the maximum value the distribution can take, in this example 20 days.

Why a gamma distrubution?

The incubation period has to be positive in value. Therefore we must specific a distribution in dist_spec which is for positive values only.

dist_spec() supports log normal and gamma distributions, which are distributions for positive values only.

For all types of delay, we will need to use distributions for positive values only - we don’t want to include delays of negative days in our analysis!

Including distribution uncertainty

To specify a variable distribution, we include uncertainty around the mean $\mu$ and standard deviation $\sigma^2$ of our gamma distribution. If our incubation period distribution has a mean $\mu$ and standard deviation $\sigma^2$, then we assume the mean ($\mu$) follows a Normal distribution with standard deviation $\sigma_{\mu}^2$:

\[\mbox{Normal}(\mu,\sigma_{\mu}^2)\]

and a standard deviation ($\sigma^2$) follows a Normal distribution with standard deviation $\sigma_{\sigma^2}^2$:

\[\mbox{Normal}(\sigma^2,\sigma_{\sigma^2}^2).\]

We specify this using dist_spec with the additional arguments mean_sd ($\sigma_{\mu}^2$) and sd_sd ($\sigma_{\sigma^2}^2$).

R

incubation_period_variable <- dist_spec(
  mean = 4, sd = 2,
  mean_sd = 0.5, sd_sd = 0.5,
  max = 20, distribution = "gamma"
)
incubation_period_variable

OUTPUT


  Uncertain gamma distribution with (untruncated) mean 4 (SD 0.5) and SD 2 (SD 0.5)

Reporting delays

After the incubation period, there will be an additional delay of time from symptom onset to case notification: the reporting delay. We can specify this as a fixed or variable distribution, or estimate a distribution from data.

When specifying a distribution, it is useful to visualise the probability density to see the peak and spread of the distribution, in this case we will use a log normal distribution. We can use the functions convert_to_logmean() and convert_to_logsd() to convert the mean and standard deviation of a normal distribution to that of a log normal distribution.

If we want to assume that the mean reporting delay is 2 days (with a standard deviation of 1 day), the log normal distribution will look like:

R

log_mean <- convert_to_logmean(2, 1)
log_sd <- convert_to_logsd(2, 1)
x <- seq(from = 0, to = 10, length = 1000)
df <- data.frame(x = x, density = dlnorm(x, meanlog = log_mean, sdlog = log_sd))
ggplot(df) +
  geom_line(
    aes(x, density)
  ) +
  theme_grey(
    base_size = 15
  )

Using the mean and standard deviation for the log normal distribution, we can specify a fixed or variable distribution using dist_spec() as before:

R

reporting_delay_variable <- dist_spec(
  mean = log_mean, sd = log_sd,
  mean_sd = 0.5, sd_sd = 0.5,
  max = 10, distribution = "lognormal"
)

If data is available on the time between symptom onset and reporting, we can use the function estimate_delay() to estimate a log normal distribution from a vector of delays. The code below illustrates how to use estimate_delay() with synthetic delay data.

R

delay_data <- rlnorm(500, log(5), 1) # synthetic delay data
reporting_delay <- estimate_delay(
  delay_data,
  samples = 1000,
  bootstraps = 10
)

Generation time

We also must specify a distribution for the generation time. Here we will use a log normal distribution with mean 3.6 and standard deviation 3.1 (Ganyani et al. 2020).

R

generation_time_variable <- dist_spec(
  mean = 3.6, sd = 3.1,
  mean_sd = 0.5, sd_sd = 0.5,
  max = 20, distribution = "lognormal"
)

Finding estimates

The function epinow() is a wrapper for the function estimate_infections() used to estimate cases by date of infection. The generation time distribution and delay distributions must be passed using the functions generation_time_opts() and delay_opts() respectively.

There are numerous other inputs that can be passed to epinow(), see EpiNow2::?epinow() for more detail. One optional input is to specify a log normal prior for the effective reproduction number $R_t$ at the start of the outbreak. We specify a mean and standard deviation as arguments of prior within rt_opts():

R

rt_log_mean <- convert_to_logmean(2, 1)
rt_log_sd <- convert_to_logsd(2, 1)
rt <- rt_opts(prior = list(mean = rt_log_mean, sd = rt_log_sd))

Bayesian inference using Stan

The Bayesian inference is performed using MCMC methods with the program Stan. There are a number of default inputs to the Stan functions including the number of chains and number of samples per chain (see ?EpiNow2::stan_opts()).

To reduce computation time, we can run chains in parallel. To do this, we must set the number of cores to be used. By default, 4 MCMC chains are run (see stan_opts()$chains), so we can set an equal number of cores to be used in parallel as follows:

R

withr::local_options(list(mc.cores = 4))

To find the maximum number of available cores on your machine, use parallel::detectCores().

Note : in the code below fixed distributions are used instead of variable. This is to speed up computation time. It is generally recommended to use variable distributions that account for additional uncertainty.

R

reported_cases <- cases[1:90, ]
estimates <- epinow(
  reported_cases = reported_cases,
  generation_time = generation_time_opts(generation_time_fixed),
  delays = delay_opts(incubation_period_fixed + reporting_delay_fixed),
  rt = rt_opts(prior = list(mean = rt_log_mean, sd = rt_log_sd))
)

OUTPUT

WARN [2024-09-24 01:23:52] epinow: There were 3 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them. -
WARN [2024-09-24 01:23:52] epinow: Examine the pairs() plot to diagnose sampling problems
 -

Results

We can extract and visualise estimates of the effective reproduction number through time:

R

estimates$plots$R

The uncertainty in the estimates increases through time. This is because estimates are informed by data in the past - within the delay periods. This difference in uncertainty is categorised into Estimate (green) utilises all data and Estimate based on partial data (orange) estimates that are based on less data (because infections that happened at the time are more likely to not have been observed yet) and therefore have increasingly wider intervals towards the date of the last data point. Finally, the Forecast (purple) is a projection ahead of time.

We can also visualise the growth rate estimate through time:

R

estimates$plots$growth_rate

To extract a summary of the key transmission metrics at the latest date in the data:

R

summary(estimates)

OUTPUT

                                 measure                 estimate
                                  <char>                   <char>
1: New confirmed cases by infection date     7267 (4048 -- 12834)
2:        Expected change in daily cases        Likely decreasing
3:            Effective reproduction no.        0.9 (0.58 -- 1.4)
4:                        Rate of growth -0.014 (-0.063 -- 0.042)
5:          Doubling/halving time (days)          -51 (16 -- -11)

As these estimates are based on partial data, they have a wide uncertainty interval.

From the summary of our analysis we see that the expected change in daily cases is Likely decreasing with the estimated new confirmed cases 7267 (4048 – 12834).
The effective reproduction number $R_t$ estimate (on the last date of the data) is 0.9 (0.58 – 1.4).
The exponential growth rate of case numbers is -0.014 (-0.063 – 0.042).
The doubling time (the time taken for case numbers to double) is -51 (16 – -11).

`Expected change in daily cases`

A factor describing expected change in daily cases based on the posterior probability that $R_t < 1$.

Probability ($p$)	Expected change
$p < 0.05$	Increasing
$0.05 \leq p< 0.4$	Likely increasing
$0.4 \leq p< 0.6$	Stable
$0.6 \leq p < 0.95$	Likely decreasing
$0.95 \leq p \leq 1$	Decreasing

Quantify geographical heterogeneity

The outbreak data of the start of the COVID-19 pandemic from the United Kingdom from the R package incidence2 includes the region in which the cases were recorded. To find regional estimates of the effective reproduction number and cases, we must format the data to have three columns:

date : the date,
region : the region,
confirm : number of confirmed cases for a region on a given date.

R

regional_cases <-
  incidence2::covidregionaldataUK[, c("date", "cases_new", "region")]
colnames(regional_cases) <- c("date", "confirm", "region")

# extract the first 90 dates for all regions
dates <- sort(unique(regional_cases$date))[1:90]
regional_cases <- regional_cases[which(regional_cases$date %in% dates), ]

head(regional_cases)

OUTPUT

        date confirm          region
1 2020-01-30      NA   East Midlands
2 2020-01-30      NA East of England
3 2020-01-30       2         England
4 2020-01-30      NA          London
5 2020-01-30      NA      North East
6 2020-01-30      NA      North West

To find regional estimates, we use the same inputs as epinow() to the function regional_epinow():

R

estimates_regional <- regional_epinow(
  reported_cases = regional_cases,
  generation_time = generation_time_opts(generation_time_fixed),
  delays = delay_opts(incubation_period_fixed + reporting_delay_fixed),
  rt = rt_opts(prior = list(mean = rt_log_mean, sd = rt_log_sd))
)

OUTPUT

INFO [2024-09-24 01:23:57] Producing following optional outputs: regions, summary, samples, plots, latest
INFO [2024-09-24 01:23:57] Reporting estimates using data up to: 2020-04-28
INFO [2024-09-24 01:23:57] No target directory specified so returning output
INFO [2024-09-24 01:23:57] Producing estimates for: East Midlands, East of England, England, London, North East, North West, Northern Ireland, Scotland, South East, South West, Wales, West Midlands, Yorkshire and The Humber
INFO [2024-09-24 01:23:57] Regions excluded: none
INFO [2024-09-24 02:07:37] Completed regional estimates
INFO [2024-09-24 02:07:37] Regions with estimates: 13
INFO [2024-09-24 02:07:37] Regions with runtime errors: 0
INFO [2024-09-24 02:07:37] Producing summary
INFO [2024-09-24 02:07:37] No summary directory specified so returning summary output
INFO [2024-09-24 02:07:38] No target directory specified so returning timings

R

estimates_regional$summary$summarised_results$table

OUTPUT

                      Region New confirmed cases by infection date
                      <char>                                <char>
 1:            East Midlands                      343 (213 -- 560)
 2:          East of England                      543 (331 -- 863)
 3:                  England                   3564 (2190 -- 5511)
 4:                   London                      293 (190 -- 452)
 5:               North East                      254 (147 -- 431)
 6:               North West                      557 (330 -- 877)
 7:         Northern Ireland                         43 (22 -- 84)
 8:                 Scotland                      283 (157 -- 539)
 9:               South East                      596 (354 -- 965)
10:               South West                      416 (293 -- 601)
11:                    Wales                        95 (64 -- 141)
12:            West Midlands                      270 (141 -- 486)
13: Yorkshire and The Humber                      481 (278 -- 784)
    Expected change in daily cases Effective reproduction no.
                            <fctr>                     <char>
 1:              Likely increasing          1.2 (0.85 -- 1.6)
 2:              Likely increasing          1.2 (0.83 -- 1.6)
 3:              Likely decreasing         0.91 (0.64 -- 1.2)
 4:              Likely decreasing         0.78 (0.56 -- 1.1)
 5:              Likely decreasing         0.91 (0.61 -- 1.3)
 6:              Likely decreasing         0.86 (0.58 -- 1.2)
 7:              Likely decreasing           0.64 (0.38 -- 1)
 8:              Likely decreasing          0.9 (0.58 -- 1.4)
 9:                         Stable         0.99 (0.67 -- 1.4)
10:                     Increasing           1.4 (1.1 -- 1.8)
11:                     Decreasing        0.57 (0.42 -- 0.77)
12:              Likely decreasing         0.71 (0.42 -- 1.1)
13:                         Stable            1 (0.69 -- 1.4)
                Rate of growth Doubling/halving time (days)
                        <char>                       <char>
 1:    0.023 (-0.021 -- 0.068)               30 (10 -- -33)
 2:    0.023 (-0.023 -- 0.067)               31 (10 -- -30)
 3:    -0.011 (-0.053 -- 0.03)              -61 (23 -- -13)
 4:   -0.03 (-0.066 -- 0.0098)              -23 (71 -- -11)
 5:   -0.011 (-0.057 -- 0.037)              -61 (19 -- -12)
 6:   -0.019 (-0.062 -- 0.023)              -37 (30 -- -11)
 7:    -0.052 (-0.1 -- 0.0053)            -13 (130 -- -6.8)
 8:   -0.013 (-0.062 -- 0.047)              -53 (15 -- -11)
 9: -0.00095 (-0.047 -- 0.046)             -730 (15 -- -15)
10:     0.046 (0.013 -- 0.084)               15 (8.3 -- 54)
11:  -0.065 (-0.094 -- -0.032)            -11 (-21 -- -7.4)
12:   -0.041 (-0.092 -- 0.012)             -17 (59 -- -7.5)
13:   0.0033 (-0.044 -- 0.051)              210 (14 -- -16)

R

estimates_regional$summary$plots$R

Summary

EpiNow2 can be used to estimate transmission metrics from case data at the start of an outbreak. The reliability of these estimates depends on the quality of the data and appropriate choice of delay distributions. In the next tutorial we will learn how to make forecasts and investigate some of the additional inference options available in EpiNow2.

Key Points

Transmission metrics can be estimated from case data after accounting for delays
Uncertainty can be accounted for in delay distributions

Content from Create a short-term forecast

Last updated on 2024-09-24 | Edit this page

Estimated time: 60 minutes

Overview

Questions

How do I create short-term forecasts from case data?
How do I account for incomplete reporting in forecasts?

Objectives

Learn how to make forecasts of cases using R package EpiNow2
Learn how to include an observation process in estimation

Prerequisites

Complete tutorial Quantifying transmission

Learners should familiarise themselves with following concept dependencies before working through this tutorial:

Statistics : probability distributions, principle of Bayesian analysis.

Epidemic theory : Effective reproduction number.

Introduction

Given case data, we can create estimates of the current and future number of cases by accounting for both delays in reporting and under reporting. To make statements about the future, we need to make an assumption of how observations up to today are related to what we expect to happen in the future. The simplest way of doing so is to assume “no change”, i.e. the reproduction number remains the same in the future as last observed. In this tutorial we will create short-term forecasts by assuming the reproduction number will remain the same as it was estimated to be on the final date for which data was available.

Create a short-term forecast

The function epinow() described in the previous tutorial is a wrapper for the function estimate_infections() used to estimate cases by date of infection. Using the same code in the previous tutorial we can extract the short-term forecast using:

R

reported_cases <- cases[1:90, ]
estimates <- epinow(
  reported_cases = reported_cases,
  generation_time = generation_time_opts(generation_time_fixed),
  delays = delay_opts(incubation_period_fixed + reporting_delay_fixed),
  rt = rt_opts(prior = list(mean = rt_log_mean, sd = rt_log_sd))
)

OUTPUT

WARN [2024-09-24 02:14:50] epinow: There were 6 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them. -
WARN [2024-09-24 02:14:50] epinow: Examine the pairs() plot to diagnose sampling problems
 -
WARN [2024-09-24 02:14:52] epinow: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
Running the chains for more iterations may help. See
https://mc-stan.org/misc/warnings.html#bulk-ess -

We can visualise the estimates of the effective reproduction number and the estimated number of cases using plot(). The estimates are split into three categories:

Estimate (green): utilises all data,
Estimate based on partial data (orange): estimates that are based less data are therefore more uncertain,
Forecast (purple): forecasts into the future.

R

plot(estimates)

Forecasting with estimates of $R_t$

By default, the short-term forecasts are created using the latest estimate of the reproduction number $R_t$. As this estimate is based on partial data, it has considerable uncertainty.

The reproduction number that is projected into the future can be changed to a less recent estimate based on more data using rt_opts():

R

rt_opts(future = "estimate")

The result will be less uncertain forecasts (as they are based on $R_t$ with a narrower uncertainty interval) but the forecasts will be based on less recent estimates of $R_t$ and assume no change since then.

Additionally, there is the option to project the value of $R_t$ into the future using a generic model by setting future = "project". As this option uses a model to forecast the value of $R_t$, the result will be forecasts that are more uncertain than estimate, for an example see here.

Incomplete observations

In the previous tutorial we accounted for delays in reporting. In EpiNow2 we also can account for incomplete observations as in reality, 100% of cases are not reported.

We will pass another input into epinow() called obs to define an observation model. The format of obs must be the obs_opt() function (see ?EpiNow2::obs_opts for more detail).

Let’s say we believe the COVID-19 outbreak data from the previous tutorial do not include all reported cases. We believe that we only observe 40% of cases. To specify this in the observation model, we must pass a scaling factor with a mean and standard deviation. If we assume that 40% of cases are in the case data (with standard deviation 1%), then we specify the scale input to obs_opts() as follows:

R

obs_scale <- list(mean = 0.4, sd = 0.01)

To run the inference framework with this observation process, we add obs = obs_opts(scale = obs_scale) to the input arguments of epinow():

R

obs_scale <- list(mean = 0.4, sd = 0.01)
reported_cases <- cases[1:90, ]
estimates <- epinow(
  reported_cases = reported_cases,
  generation_time = generation_time_opts(generation_time_fixed),
  delays = delay_opts(incubation_period_fixed + reporting_delay_fixed),
  rt = rt_opts(prior = list(mean = rt_log_mean, sd = rt_log_sd)),
  obs = obs_opts(scale = obs_scale)
)

OUTPUT

WARN [2024-09-24 02:21:47] epinow: There were 4 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them. -
WARN [2024-09-24 02:21:47] epinow: Examine the pairs() plot to diagnose sampling problems
 -

R

summary(estimates)

OUTPUT

                                 measure                 estimate
                                  <char>                   <char>
1: New confirmed cases by infection date   18060 (10134 -- 32193)
2:        Expected change in daily cases        Likely decreasing
3:            Effective reproduction no.        0.9 (0.57 -- 1.3)
4:                        Rate of growth -0.014 (-0.065 -- 0.039)
5:          Doubling/halving time (days)          -51 (18 -- -11)

The estimates of transmission measures such as the effective reproduction number and rate of growth are similar (or the same in value) compared to when we didn’t account for incomplete observations (see previous tutorial). However the number of new confirmed cases by infection date has changed substantially in magnitude to reflect the assumption that only 40% of cases are in the data set.

We can also change the default distribution from Negative Binomial to Poisson, remove the default week effect and more. See ?EpiNow2::obs_opts for more detail.

Forecast secondary observations

EpiNow2 also has the ability to estimate and forecast secondary observations e.g. deaths, hospitalisations from a primary observation e.g. cases. Here we will illustrate how to forecast the number of deaths arising from observed cases of COVID-19 in the early stages of the UK outbreak.

First, we must format our data to have the following columns:

date : the date (as a date object see ?is.Date()),
primary : number of primary observations on that date, in this example cases,
secondary : number of secondary observations date, in this example deaths.

R

reported_cases_deaths <- aggregate(
  cbind(cases_new, deaths_new) ~ date,
  data =
    incidence2::covidregionaldataUK[, c("date", "cases_new", "deaths_new")],
  FUN = sum
)
colnames(reported_cases_deaths) <- c("date", "primary", "secondary")

Using the first 30 days of data on cases and deaths, we will estimate the relationship between the primary and secondary observations using estimate_secondary(), then forecast future deaths using forecast_secondary(). For detail on the model see the model documentation.

We must specify the type of observation using type in secondary_opts(), options include:

“incidence” : secondary observations arise from previous primary observations, i.e. deaths arising from recorded cases.
“prevalence” : secondary observations arise from a combination current primary observations and past secondary observations, i.e. hospital bed usage arising from current hospital admissions and past hospital bed usage.

In this example we specify secondary_opts(type = "incidence"). See ?EpiNow2::secondary_opts for more detail).

The final key input is the delay distribution between the primary and secondary observations. Here this is the delay between case report and death, we assume this follows a gamma distribution with mean of 14 days and standard deviation of 5 days. Using dist_spec() we specify a fixed gamma distribution.

There are further function inputs to estimate_secondary() which can be specified, including adding an observation process, see ?EpiNow2::estimate_secondary for detail on the options.

To find the model fit between cases and deaths :

R

estimate_cases_to_deaths <- estimate_secondary(
  reports = reported_cases_deaths[1:30, ],
  secondary = secondary_opts(type = "incidence"),
  delays = delay_opts(dist_spec(
    mean = 14, sd = 5,
    max = 30, distribution = "gamma"
  ))
)

Be cautious of timescale

In the early stages of an outbreak there can be substantial changes in testing and reporting. If there are testing changes from one month to another, then there will be a bias in the model fit. Therefore, you should be cautious of the time-scale of data used in the model fit and forecast.

We plot the model fit (shaded ribbons) with the secondary observations (bar plot) and primary observations (dotted line) as follows:

R

plot(estimate_cases_to_deaths, primary = TRUE)

To use this model fit to forecast deaths, we pass a data frame consisting of the primary observation (cases) for dates not used in the model fit.

Note : in this tutorial we are using data where we know the deaths and cases, so we create a data frame by extracting the cases. But in practice, this would be a different data set consisting of cases only.

R

cases_to_forecast <- reported_cases_deaths[31:60, c("date", "primary")]
colnames(cases_to_forecast) <- c("date", "value")

To forecast, we use the model fit estimate_cases_to_deaths:

R

deaths_forecast <- forecast_secondary(
  estimate = estimate_cases_to_deaths,
  primary = cases_to_forecast
)
plot(deaths_forecast)

WARNING

Warning: Removed 30 rows containing missing values (`position_stack()`).

The plot shows the forecast secondary observations (deaths) over the dates which we have recorded cases for. It is also possible to forecast deaths using forecast cases, here you would specify primary as the estimates output from estimate_infections().

Challenge : Ebola outbreak analysis

Challenge

Download the file ebola_cases.csv and read it into R. The simulated data consists of the date of symptom onset and number of confirmed cases of the early stages of the Ebola outbreak in Sierra Leone in 2014.

Using the first 3 months (120 days) of data:

Estimate of cases increasing or decreasing on day 120 of the outbreak (Hint: Find the effective reproduction number and growth rate on day 120)
Create a two week forecast of number of cases

You can use the following parameter values for the delay distribution(s) and generation time distribution.

Incubation period : Log normal$(2.487,0.330)$ (Eichner et al. 2011 via {epiparameter})
Generation time : Gamma$(15.3, 10.1)$ (WHO Ebola Response Team 2014)

You may include some uncertainty around the mean and standard deviation of these distributions.

HINT : data format

Ensure the data is in the correct format :

date : the date (as a date object see ?is.Date()),
confirm : number of confirmed cases on that date.

SOLUTION

To estimate the effective reproduction number and growth rate, we will use the function epinow().

As the data consists of date of symptom onset, we only need to specify a delay distribution for the incubation period and the generation time.

We specify the distributions with some uncertainty around the mean and standard deviation of the log normal distribution for the incubation period and the Gamma distribution for the generation time.

R

ebola_incubation_period <- dist_spec(
  mean = 2.487, sd = 0.330,
  mean_sd = 0.5, sd_sd = 0.5,
  max = 20, distribution = "lognormal"
)

ebola_generation_time <- dist_spec(
  mean = 15.3, sd = 10.1,
  mean_sd = 0.5, sd_sd = 0.5,
  max = 30, distribution = "gamma"
)

As we want to also create a two week forecast, we specify horizon = 14 to forecast 14 days instead of the default 7 days.

R

# read data
# e.g.: if path to file is data/raw-data/ebola_cases.csv then:
ebola_cases <-
  read.csv(here::here("data", "raw-data", "ebola_cases.csv"))

R

# format date column
ebola_cases$date <- as.Date(ebola_cases$date)

ebola_estimates <- epinow(
  reported_cases = ebola_cases[1:120, ], # first 3 months of data only
  generation_time = generation_time_opts(ebola_generation_time),
  delays = delay_opts(ebola_incubation_period),
  # horizon needs to be 14 days to create two week forecast (default is 7 days)
  horizon = 14
)

OUTPUT

WARN [2024-09-24 02:23:54] epinow: There were 24 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them. -
WARN [2024-09-24 02:23:54] epinow: Examine the pairs() plot to diagnose sampling problems
 -

R

summary(ebola_estimates)

OUTPUT

                                 measure                estimate
                                  <char>                  <char>
1: New confirmed cases by infection date         102 (47 -- 265)
2:        Expected change in daily cases              Increasing
3:            Effective reproduction no.        1.7 (1.1 -- 2.9)
4:                        Rate of growth 0.043 (0.0041 -- 0.091)
5:          Doubling/halving time (days)         16 (7.6 -- 170)

The effective reproduction number $R_t$ estimate (on the last date of the data) is 1.7 (1.1 – 2.9). The exponential growth rate of case numbers is 0.043 (0.0041 – 0.091).

Visualize the estimates:

R

plot(ebola_estimates)

Summary

EpiNow2 can be used to create short term forecasts and to estimate the relationship between different outcomes. There are a range of model options that can be implemented for different analysis, including adding an observational process to account for incomplete reporting. See the vignette for more details on different model options in EpiNow2 that aren’t covered in these tutorials.

Key Points

We can create short-term forecasts by making assumptions about the future behaviour of the reproduction number
Incomplete case reporting can be accounted for in estimates

Content from Simulating transmission

Last updated on 2024-09-24 | Edit this page

Estimated time: 75 minutes

Overview

Questions

How do I simulate disease spread using a mathematical model?
What inputs are needed for a model simulation?
How do I account for uncertainty?

Objectives

Load an existing model structure from {epidemics} R package
Load an existing social contact matrix with socialmixr
Generate a disease spread model simulation with {epidemics}
Generate multiple model simulations and visualise uncertainty

Prerequisites

Learners should familiarise themselves with following concept dependencies before working through this tutorial:

Mathematical Modelling : Introduction to infectious disease models, state variables, model parameters, initial conditions, ordinary differential equations.

Epidemic theory : Transmission, Reproduction number.

Introduction

Mathematical models are useful tools for generating future trajectories of disease spread. In this tutorial, we will use the R package {epidemics} to generate disease trajectories of an influenza strain with pandemic potential. By the end of this tutorial, you will be able to generate the trajectory below showing the number of infectious individuals in different age categories over time.

Instructor Note

By the end of this tutorial, learners should be able to replicate the above image on their own computers.

Simulating disease spread

To generate predictions of infectious disease trajectories, we must first select a mathematical model to use. There is a library of models to choose from in epidemics. Models in epidemics are prefixed with model and suffixed by the name of infection (e.g. Ebola) or a different identifier (e.g. default), and whether the model has a R or C++ code base.

In this tutorial, we will use the default model in epidemics, model_default_cpp() which is an age-structured SEIR model described by a system of ordinary differential equations. For each age group $i$, individuals are classed as either susceptible $S$, infected but not yet infectious $E$, infectious $I$ or recovered $R$. The schematic below shows the processes which describe the flow of individuals between the disease states $S$, $E$, $I$ and $R$ and the key parameters for each process.

Model parameters : rates

In ODE models, model parameters are often (but not always) specified as rates. The rate at which an event occurs is the inverse of the average time until that event. For example, in the SEIR model, the recovery rate $\gamma$ is the inverse of the average infectious period.

We can use knowledge of the natural history of the disease to inform our values of rates. If the average infectious period of an infection is 8 days, then the daily recovery rate is $\gamma = 1/8 = 0.125$.

For each disease state ($S$, $E$, $I$ and $R$) and age group ($i$), we have an ordinary differential equation describing the rate of change with respect to time.

\[ \begin{aligned} \frac{dS_i}{dt} & = - \beta S_i \sum_j C_{i,j} I_j \\ \frac{dE_i}{dt} &= \beta S_i\sum_j C_{i,j} I_j - \alpha E_i \\ \frac{dI_i}{dt} &= \alpha E_i - \gamma I_i \\ \frac{dR_i}{dt} &=\gamma I_i \\ \end{aligned} \] Individuals in age group ($i$) move from the susceptible state ($S_i$) to the exposed state ($E_i$) via age group specific contact with the infectious individuals in their own and other age groups $\beta S_i \sum_j C_{i,j} I_j$. The contact matrix $C$ allows for heterogeneity in contacts between age groups. They then move to the infectious state at a rate $\alpha$ and recover at a rate $\gamma$. There is no loss of immunity (there are no flows out of the recovered state).

The model parameters definitions are :

transmission rate or transmissibility $\beta$,
contact matrix $C$ containing the frequency of contacts between age groups (a square $i \times j$ matrix),
infectiousness rate $\alpha$ (preinfectious period (latent period) =$1/\alpha$),
recovery rate $\gamma$ (infectious period = $1/\gamma$).

Exposed, infected, infectious

Confusion sometimes arises when referring to the terms ‘exposed’, ‘infected’ and ‘infectious’ in mathematical modelling. Infection occurs after a person has been exposed, but in modelling terms individuals that are ‘exposed’ are treated as already infected.

We will use the following definitions for our state variables:

$E$ = Exposed : infected but not yet infectious,
$I$ = Infectious: infected and infectious.

To generate trajectories using our model, we must prepare the following inputs :

Contact matrix
Initial conditions
Population structure
Model parameters

1. Contact matrix

Contact matrices can be estimated from surveys or contact data, or synthetic ones can be used. We will use the R package socialmixr to load in a contact matrix estimated from POLYMOD survey data (Mossong et al. 2008).

Load contact and population data

Using the R package socialmixr, run the following lines of R code to obtain the contact matrix for the United Kingdom for the year age bins:

age between 0 and 20 years,
age between 20 and 40,
40 years and over.

R

polymod <- socialmixr::polymod
contact_data <- socialmixr::contact_matrix(
  survey = polymod,
  countries = "United Kingdom",
  age.limits = c(0, 20, 40),
  symmetric = TRUE
)
# prepare contact matrix
contact_matrix <- t(contact_data$matrix)
contact_matrix

Output

OUTPUT


contact.age.group     [,1]     [,2]     [,3]
          [0,20)  7.883663 2.794154 1.565665
          [20,40) 3.120220 4.854839 2.624868
          40+     3.063895 4.599893 5.005571

The result is a square matrix with rows and columns for each age group. Contact matrices can be loaded from other sources, but they must be formatted as a matrix to be used in epidemics.

Why would a contact matrix be non-symmetric?

One of the arguments of the function contact_matrix() is symmetric=TRUE. This means that the total number of contacts of age group 1 with age group 2, should be the same as the total number of contacts of age group 2 and age group 1 (see the socialmixr vignette for more detail). However, when contact matrices are estimated from surveys or other sources, the reported number of contacts may differ by age group resulting in a non-symmetric contact matrix (Prem et al 2021).

2. Initial conditions

The initial conditions are the proportion of individuals in each disease state $S$, $E$, $I$ and $R$ for each age group at time 0. In this example, we have three age groups age between 0 and 20 years, age between 20 and 40 years and over. Let’s assume that in the youngest age category, one in a million individuals are infectious, and the remaining age categories are infection free.

The initial conditions in the first age category are $S(0)=1-\frac{1}{1,000,000}$, $E(0) =0$, $I(0)=\frac{1}{1,000,000}$, $R(0)=0$. This is specified as a vector as follows:

R

initial_i <- 1e-6
initial_conditions_inf <- c(
  S = 1 - initial_i, E = 0, I = initial_i, R = 0, V = 0
)

For the age categories that are free from infection, the initial conditions are $S(0)=1$, $E(0) =0$, $I(0)=0$, $R(0)=0$. We specify this as follows,

R

initial_conditions_free <- c(
  S = 1, E = 0, I = 0, R = 0, V = 0
)

We combine the three initial conditions vectors into one matrix,

R

# combine the initial conditions
initial_conditions <- rbind(
  initial_conditions_inf, # age group 1
  initial_conditions_free, # age group 2
  initial_conditions_free # age group 3
)

# use contact matrix to assign age group names
rownames(initial_conditions) <- rownames(contact_matrix)
initial_conditions

OUTPUT

               S E     I R V
[0,20)  0.999999 0 1e-06 0 0
[20,40) 1.000000 0 0e+00 0 0
40+     1.000000 0 0e+00 0 0

3. Population structure

The population object requires a vector containing the demographic structure of the population. The demographic vector must be a named vector containing the number of individuals in each age group of our given population. In this example, we can extract the demographic information from the contact_data object that we obtained using the socialmixr package.

R

demography_vector <- contact_data$demography$population
names(demography_vector) <- rownames(contact_matrix)
demography_vector

OUTPUT

  [0,20)  [20,40)      40+
14799290 16526302 28961159

To create our population object, we call the function population() specifying a name, the contact matrix, the demography vector and the initial conditions.

R

uk_population <- population(
  name = "UK",
  contact_matrix = contact_matrix,
  demography_vector = demography_vector,
  initial_conditions = initial_conditions
)

4. Model parameters

To run our model we need to specify the model parameters:

transmissibility $\beta$,
infectiousness rate $\alpha$ (preinfectious period=$1/\alpha$),
recovery rate $\gamma$ (infectious period=$1/\gamma$).

In epidemics, we specify the model inputs as :

transmissibility = $R_0 \gamma$,
infectiousness_rate = $\alpha$,
recovery_rate = $\gamma$,

We will simulate a strain of influenza with pandemic potential with $R_0=1.46$, a preinfectious period of 3 days and infectious period of 7 days. Therefore our inputs will be:

transmissibility = 1.46 / 7.0,
infectiousness_rate = 1.0 / 3.0,
recovery_rate = 1.0 / 7.0.

The basic reproduction number $R_0$

The basic reproduction number, $R_0$, for the SEIR model is:

\[ R_0 = \frac{\beta}{\gamma}.\]

Therefore, we can rewrite transmissibility $\beta$, as:

\[ \beta = R_0 \gamma.\]

Running the model

Running (solving) the model

For models that are described by ODEs, running the model actually means to solve the system of ODEs. ODEs describe the rate of change in the disease states with respect to time, to find the number of individuals in each of these states, we use numerical integration methods to solve the equations.

In epidemics, the ODE solver uses the Runge-Kutta method.

Now we are ready to run our model. Let’s load the epidemics package :

R

library(epidemics)

Then we specify time_end=600 to run the model for 600 days.

R

output <- model_default_cpp(
  population = uk_population,
  transmissibility = 1.46 / 7.0,
  infectiousness_rate = 1.0 / 3.0,
  recovery_rate = 1.0 / 7.0,
  time_end = 600, increment = 1.0
)
head(output)

OUTPUT

  time demography_group compartment    value
1    0           [0,20) susceptible 14799275
2    0          [20,40) susceptible 16526302
3    0              40+ susceptible 28961159
4    0           [0,20)     exposed        0
5    0          [20,40)     exposed        0
6    0              40+     exposed        0

Note : This model also has the functionality to include vaccination and tracks the number of vaccinated individuals through time. Even though we have not specified any vaccination, there is still a vaccinated compartment in the output (containing no individuals). We will cover the use of vaccination in future tutorials.

Our model output consists of the number of individuals in each compartment in each age group through time. We can visualise the infectious individuals only (those in the $I$ class) through time.

R

filter(output_plot, compartment %in% c("exposed", "infectious")) %>%
  ggplot(
    aes(
      x = time,
      y = value,
      col = demography_group,
      linetype = compartment
    )
  ) +
  geom_line(
    linewidth = 1.2
  ) +
  scale_y_continuous(
    labels = scales::comma
  ) +
  scale_colour_brewer(
    palette = "Dark2",
    name = "Age group"
  ) +
  expand_limits(
    y = c(0, 500e3)
  ) +
  coord_cartesian(
    expand = FALSE
  ) +
  theme_bw(
    base_size = 15
  ) +
  theme(
    legend.position = "top"
  ) +
  labs(
    x = "Simulation time (days)",
    linetype = "Compartment",
    y = "Individuals"
  )

Time increments

Note that there is a default argument of increment = 1. This relates to the time step of the ODE solver. When the parameters are specified on a daily time-scale and maximum number of time steps (time_end) is days, the default time step of the ODE solver one day.

The choice of increment will depend on the time scale of the parameters, and the rate at which events can occur. In general, the increment should smaller than the fastest event. For example, if parameters are on a monthly time scale, but some events will occur within a month, then the increment should be less than one month.

Accounting for uncertainty

As the epidemic model is deterministic, we have one trajectory for our given parameter values. In practice, we have uncertainty in the value of our parameters. To account for this, we must run our model for different parameter combinations.

We ran our model with $R_0= 1.5$. However, we believe that $R_0$ follows a normal distribution with mean 1.5 and standard deviation 0.05. To account for uncertainty we will run the model for different values of $R_0$. The steps we will follow to do this are:

Obtain 100 samples from the from a normal distribution

R

R0_vec <- rnorm(100, 1.5, 0.05)

Run the model 100 times with $R_0$ equal to a different sample each time

R

output_samples <- Map(
  R0_vec,
  seq_along(R0_vec),
  f = function(x, i) {
    # run an epidemic model using `epidemic()`
    output <- model_default_cpp(
      population = uk_population,
      transmissibility = x / 7.0,
      infectiousness_rate = 1.0 / 3.0,
      recovery_rate = 1.0 / 7.0,
      time_end = 600, increment = 1.0
    )

    # add replicate number and return data
    output$replicate <- x
    output
  }
)

# combine to prepare for plotting
output_samples <- bind_rows(output_samples)

Calculate the mean and 95% quantiles of number of infectious individuals across each model simulation and visualise output

R

ggplot(
  output_samples[output_samples$compartment == "infectious", ],
  aes(time, value)
) +
  stat_summary(geom = "line", fun = mean) +
  stat_summary(
    geom = "ribbon",
    fun.min = function(z) {
      quantile(z, 0.025)
    },
    fun.max = function(z) {
      quantile(z, 0.975)
    },
    alpha = 0.3
  ) +
  facet_grid(
    cols = vars(demography_group)
  ) +
  labs(
    x = "Simulation time (days)",
    y = "Individuals"
  ) +
  theme_bw(
    base_size = 15
  )

Deciding which parameters to include uncertainty in depends on a few factors: how well informed a parameter value is e.g. consistency of estimates from the literature; how sensitive model outputs are to parameter value changes; and the purpose of the modelling task. See McCabe et al. 2021 to learn about different types of uncertainty in infectious disease modelling.

Summary

In this tutorial, we have learnt how to simulate disease spread using a mathematical model. Once a model has been chosen, the parameters and other inputs must be specified in the correct way to perform model simulations. In the next tutorial, we will consider how to choose the right model for different tasks.

Key Points

Disease trajectories can be generated using the R package epidemics
Uncertainty should be included in model trajectories using a range of model parameter values

Content from Choosing an appropriate model

Last updated on 2024-09-24 | Edit this page

Estimated time: 30 minutes

Overview

Questions

How do I choose a mathematical model that’s appropriate to complete my analytical task?

Objectives

Understand the model requirements for a specific research question

Prerequisites

Complete tutorial Simulating transmission

Introduction

There are existing mathematical models for different infections, interventions and transmission patterns which can be used to answer new questions. In this tutorial, we will learn how to choose an existing model to complete a given task.

Instructor Note

The focus of this tutorial is understanding existing models to decide if they are appropriate for a defined question.

Choosing a model

When deciding which mathematical model to use, there are a number of questions we must consider :

What is the infection/disease of interest?

A model may already exist for your study disease, or there may be a model for an infection that has the same transmission pathways and epidemiological features that can be used.

Do we need a deterministic or stochastic model?

Model structures differ for whether the disease has pandemic potential or not. When predicted numbers of infection are small, stochastic variation in output can have an effect on whether an outbreak takes off or not. Outbreaks are usually smaller in magnitude than epidemics, so its often appropriate to use a stochastic model to characterise the uncertainty in the early stages of the outbreak. Epidemics are larger in magnitude than outbreaks and so a deterministic model is suitable as we have less interest in the stochastic variation in output.

What is the outcome of interest?

The outcome of interest can be a feature of a mathematical model. It may be that you are interested in the predicted numbers of infection through time, or in a specific outcome such as hospitalisations or cases of severe disease.

How is transmission modelled?

For example, direct or indirect, airborne or vector-borne.

How are the different processes (e.g. transmission) formulated in the equations?

There can be subtle differences in model structures for the same infection or outbreak type which can be missed without studying the equations. For example, transmissibility parameters can be specified as rates or probabilities. If you want to use parameter values from other published models, you must check that transmission is formulated in the same way.

Will any interventions be modelled?

Finally, interventions such as vaccination may be of interest. A model may or may not have the capability to include the impact of different interventions on different time scales (continuous time or at discrete time points). We discuss interventions in detail in the tutorial Modelling interventions.

Available models in `epidemics`

The R package epidemics contains functions to run existing models. For details on the models that are available, see the package vignettes. To learn how to run the models in R, read the documentation using ?epidemics::model_ebola_r.

What model?

You have been asked to explore the variation in numbers of infectious individuals in the early stages of an Ebola outbreak.

Which of the following models would be an appropriate choice for this task:

model_default_cpp()
model_ebola_r()

HINT

Consider the following questions:

Checklist

What is the infection/disease of interest?
Do we need a deterministic or stochastic model?
What is the outcome of interest?
Will any interventions be modelled?

SOLUTION

What is the infection/disease of interest? Ebola
Do we need a deterministic or stochastic model? A stochastic model would allow us to explore variation in the early stages of the outbreak
What is the outcome of interest? Number of infections
Will any interventions be modelled? No

`model_default_cpp()`

A deterministic SEIR model with age specific direct transmission.

The model is capable of predicting an Ebola type outbreak, but as the model is deterministic, we are not able to explore stochastic variation in the early stages of the outbreak.

`model_ebola_r()`

A stochastic SEIHFR (Susceptible, Exposed, Infectious, Hospitalised, Funeral, Removed) model that was developed specifically for infection with Ebola. The model has stochasticity in the passage times between states, which are modelled as Erlang distributions.

The key parameters affecting the transition between states are:

$R_0$, the basic reproduction number,
$\rho^I$, the mean infectious period,
$\rho^E$, the mean preinfectious period,
$p_{hosp}$ the probability of being transferred to the hospitalised compartment.

Note: the functional relationship between the preinfectious period ($\rho^E$) and the transition rate between exposed and infectious ($\gamma^E$) is $\rho^E = k^E/\gamma^E$ where $k^E$ is the shape of the Erlang distribution. Similarly for the infectious period $\rho^I = k^I/\gamma^I$. See here for more detail on the stochastic model formulation.

The model has additional parameters describing the transmission risk in hospital and funeral settings:

$p_{ETU}$, the proportion of hospitalised cases contributing to the infection of susceptibles (ETU = Ebola virus treatment units),
$p_{funeral}$, the proportion of funerals at which the risk of transmission is the same as of infectious individuals in the community.

As this model is stochastic, it is the most appropriate choice to explore how variation in numbers of infected individuals in the early stages of an Ebola outbreak.

Challenge : Ebola outbreak analysis

Running the model

You have been tasked to generate initial trajectories of an Ebola outbreak in Guinea. Using model_ebola_r() and the the information detailed below, complete the following tasks:

Run the model once and plot the number of infectious individuals through time
Run model 100 times and plot the mean, upper and lower 95% quantiles of the number of infectious individuals through time

Population size : 14 million
Initial number of exposed individuals : 10
Initial number of infectious individuals : 5
Time of simulation : 120 days
Parameter values :
- $R_0$ (r0) = 1.1,
- $p^I$ (infectious_period) = 12,
- $p^E$ (preinfectious_period) = 5,
- $k^E=k^I = 2$,
- $1-p_{hosp}$ (prop_community) = 0.9,
- $p_{ETU}$ (etu_risk) = 0.7,
- $p_{funeral}$ (funeral_risk) = 0.5

Code for initial conditions

R

# set population size
population_size <- 14e6

E0 <- 10
I0 <- 5
# prepare initial conditions as proportions
initial_conditions <- c(
  S = population_size - (E0 + I0), E = E0, I = I0, H = 0, F = 0, R = 0
) / population_size

guinea_population <- population(
  name = "Guinea",
  contact_matrix = matrix(1), # note dummy value
  demography_vector = population_size, # 14 million, no age groups
  initial_conditions = matrix(
    initial_conditions,
    nrow = 1
  )
)

HINT : Multiple model simulations

Adapt the code from the accounting for uncertainty section

SOLUTION

Run the model once and plot the number of infectious individuals through time

R

output <- model_ebola_r(
  population = guinea_population,
  transmissibility = 1.1 / 12,
  infectiousness_rate = 2.0 / 5,
  removal_rate = 2.0 / 12,
  prop_community = 0.9,
  etu_risk = 0.7,
  funeral_risk = 0.5,
  time_end = 100
)

ggplot(output[output$compartment == "infectious", ]) +
  geom_line(
    aes(time, value),
    linewidth = 1.2
  ) +
  scale_y_continuous(
    labels = scales::comma
  ) +
  labs(
    x = "Simulation time (days)",
    y = "Individuals"
  ) +
  theme_bw(
    base_size = 15
  )

Run model 100 times and plot the mean, upper and lower 95% quantiles of the number of infectious individuals through time

We run the model 100 times with the same parameter values.

R

output_samples <- Map(
  x = seq(100),
  f = function(x) {
    output <- model_ebola_r(
      population = guinea_population,
      transmissibility = 1.1 / 12,
      infectiousness_rate = 2.0 / 5,
      removal_rate = 2.0 / 12,
      prop_community = 0.9,
      etu_risk = 0.7,
      funeral_risk = 0.5,
      time_end = 100
    )
    # add replicate number and return data
    output$replicate <- x
    output
  }
)

output_samples <- bind_rows(output_samples) # requires the dplyr package

ggplot(
  output_samples[output_samples$compartment == "infectious", ],
  aes(time, value)
) +
  stat_summary(geom = "line", fun = mean) +
  stat_summary(
    geom = "ribbon",
    fun.min = function(z) {
      quantile(z, 0.025)
    },
    fun.max = function(z) {
      quantile(z, 0.975)
    },
    alpha = 0.3
  ) +
  labs(
    x = "Simulation time (days)",
    y = "Individuals"
  ) +
  theme_bw(
    base_size = 15
  )

Key Points

Existing mathematical models should be selected according to the research question
It is important to check that a model has appropriate assumptions about transmission, outbreak potential, outcomes and interventions

Content from Modelling interventions

Last updated on 2024-09-24 | Edit this page

Estimated time: 75 minutes

Overview

Questions

How do I investigate the effect of interventions on disease trajectories?

Objectives

Add pharmaceutical and non-pharmaceutical interventions to an {epidemics} model

Prerequisites

Complete tutorial Simulating transmission

Learners should familiarise themselves with following concept dependencies before working through this tutorial:

Outbreak response : Intervention types.

Introduction

Mathematical models can be used to generate trajectories of disease spread under the implementation of interventions at different stages of an outbreak. These predictions can be used to make decisions on what interventions could be implemented to slow down the spread of diseases.

We can assume interventions in mathematical models reduce the values of relevant parameters e.g. reduce transmissibility while in place. Or it may be appropriate to assume individuals are classified into a new disease state, e.g. once vaccinated we assume individuals are no longer susceptible to infection and therefore move to a vaccinated state. In this tutorial, we will introduce how to include three different interventions in model of COVID-19 transmission.

Instructor Note

In this tutorial different types of intervention and how they can be modelled are introduced. Learners should be able to understand the underlying mechanism of these interventions (e.g. reduce contact rate) as well as how to implement the code to include such interventions.

Non-pharmaceutical interventions

Non-pharmaceutical interventions (NPIs) are measures put in place to reduce transmission that do not include the administration of drugs or vaccinations. NPIs aim reduce contact between infectious and susceptible individuals. For example, washing hands, wearing masks and closures of school and workplaces.

We will investigate the effect of interventions on a COVID-19 outbreak using an SEIR model (model_default_cpp() in the R package {epidemics}). We will set $R_0 = 2.7$, latent period or preinfectious period $= 4$ and the infectious_period $= 5.5$ (parameters adapted from Davies et al. (2020)). We load a contact matrix with age bins 0-18, 18-65, 65 years and older using socialmixr and assume that one in every 1 million in each age group is infectious at the start of the epidemic.

R

polymod <- socialmixr::polymod
contact_data <- socialmixr::contact_matrix(
  polymod,
  countries = "United Kingdom",
  age.limits = c(0, 15, 65),
  symmetric = TRUE
)

# prepare contact matrix
contact_matrix <- t(contact_data$matrix)

# prepare the demography vector
demography_vector <- contact_data$demography$population
names(demography_vector) <- rownames(contact_matrix)

# initial conditions: one in every 1 million is infected
initial_i <- 1e-6
initial_conditions <- c(
  S = 1 - initial_i, E = 0, I = initial_i, R = 0, V = 0
)

# build for all age groups
initial_conditions <- matrix(
  rep(initial_conditions, dim(contact_matrix)[1]),
  ncol = 5, byrow = TRUE
)
rownames(initial_conditions) <- rownames(contact_matrix)

# prepare the population to model as affected by the epidemic
uk_population <- population(
  name = "UK",
  contact_matrix = contact_matrix,
  demography_vector = demography_vector,
  initial_conditions = initial_conditions
)

Effect of school closures on COVID-19 spread

The first NPI we will consider is the effect of school closures on reducing the number of individuals infectious with COVID-19 through time. We assume that a school closure will reduce the frequency of contacts within and between different age groups. We assume that school closures will reduce the contacts between school aged children (aged 0-15) by 0.5, and will cause a small reduction (0.01) in the contacts between adults (aged 15 and over).

To include an intervention in our model we must create an intervention object. The inputs are the name of the intervention (name), the type of intervention (contacts or rate), the start time (time_begin), the end time (time_end) and the reduction (reduction). The values of the reduction matrix are specified in the same order as the age groups in the contact matrix.

R

rownames(contact_matrix)

OUTPUT

[1] "[0,15)"  "[15,65)" "65+"

Therefore, we specify reduction = matrix(c(0.5, 0.01, 0.01)). We assume that the school closures start on day 50 and are in place for a further 100 days. Therefore our intervention object is :

R

close_schools <- intervention(
  name = "School closure",
  type = "contacts",
  time_begin = 50,
  time_end = 50 + 100,
  reduction = matrix(c(0.5, 0.01, 0.01))
)

Effect of interventions on contacts

In epidemics, the contact matrix is scaled down by proportions for the period in which the intervention is in place. To understand how the reduction is calculated within the model functions, consider a contact matrix for two age groups with equal number of contacts:

OUTPUT

     [,1] [,2]
[1,]    1    1
[2,]    1    1

If the reduction is 50% in group 1 and 10% in group 2, the contact matrix during the intervention will be:

OUTPUT

     [,1] [,2]
[1,] 0.25 0.45
[2,] 0.45 0.81

The contacts within group 1 are reduced by 50% twice to accommodate for a 50% reduction in outgoing and incoming contacts ($1\times 0.5 \times 0.5 = 0.25$). Similarly, the contacts within group 2 are reduced by 10% twice. The contacts between group 1 and group 2 are reduced by 50% and then by 10% ($1 \times 0.5 \times 0.9= 0.45$).

Using transmissibility $= 2.7/5.5$ (remember that transmissibility = $R_0$/ infectious period), infectiousness rate $1/= 4$ and the recovery rate $= 1/5.5$ we run the model withintervention = list(contacts = close_schools) as follows :

R

output_school <- model_default_cpp(
  population = uk_population,
  transmissibility = 2.7 / 5.5,
  infectiousness_rate = 1.0 / 4.0,
  recovery_rate = 1.0 / 5.5,
  intervention = list(contacts = close_schools),
  time_end = 300, increment = 1.0
)

To be able to see the effect of our intervention, we also run the model where there is no intervention, combine the two outputs into one data frame and then plot the output. Here we plot the total number of infectious individuals in all age groups using ggplot2::stat_summary():

R

# run baseline simulation with no intervention
output_baseline <- model_default_cpp(
  population = uk_population,
  transmissibility = 2.7 / 5.5,
  infectiousness_rate = 1.0 / 4.0,
  recovery_rate = 1.0 / 5.5,
  time_end = 300, increment = 1.0
)

# create intervention_type column for plotting
output_school$intervention_type <- "school closure"
output_baseline$intervention_type <- "baseline"
output <- rbind(output_school, output_baseline)

ggplot(data = output[output$compartment == "infectious", ]) +
  aes(
    x = time,
    y = value,
    color = intervention_type,
    linetype = intervention_type
  ) +
  stat_summary(
    fun = "sum",
    geom = "line",
    linewidth = 1
  ) +
  scale_y_continuous(
    labels = scales::comma
  ) +
  labs(
    x = "Simulation time (days)",
    y = "Individuals"
  ) +
  theme_bw(
    base_size = 15
  ) +
  geom_vline(
    xintercept = c(close_schools$time_begin, close_schools$time_end),
    colour = "black",
    linetype = "dashed",
    linewidth = 0.2
  ) +
  annotate(
    geom = "text",
    label = "Schools closed",
    colour = "black",
    x = (close_schools$time_end - close_schools$time_begin) / 2 +
      close_schools$time_begin,
    y = 10,
    angle = 0,
    vjust = "outward"
  )

We see that with the intervention in place, the infection still spreads through the population, though the peak number of infectious individuals is smaller than the baseline with no intervention in place (solid line).

Effect of mask wearing on COVID-19 spread

We can model the effect of other NPIs as reducing the value of relevant parameters. For example, we want to investigate the effect of mask wearing on the number of individuals infectious with COVID-19 through time.

We expect that mask wearing will reduce an individual’s infectiousness. As we are using a population based model, we cannot make changes to individual behaviour and so assume that the transmissibility $\beta$ is reduced by a proportion due to mask wearing in the population. We specify this proportion, $\theta$ as product of the proportion wearing masks multiplied by the proportion reduction in transmissibility (adapted from Li et al. 2020)

We create an intervention object with type = rate and reduction = 0.161. Using parameters adapted from Li et al. 2020 we have proportion wearing masks = coverage $\times$ availability = $0.54 \times 0.525 = 0.2835$, proportion reduction in transmissibility = $0.575$. Therefore, $\theta = 0.2835 \times 0.575 = 0.163$. We assume that the mask wearing mandate starts at day 40 and is in place for 200 days.

R

mask_mandate <- intervention(
  name = "mask mandate",
  type = "rate",
  time_begin = 40,
  time_end = 40 + 200,
  reduction = 0.163
)

To implement this intervention on the parameter $\beta$, we specify intervention = list(beta = mask_mandate).

R

output_masks <- model_default_cpp(
  population = uk_population,
  transmissibility = 2.7 / 5.5,
  infectiousness_rate = 1.0 / 4.0,
  recovery_rate = 1.0 / 5.5,
  intervention = list(transmissibility = mask_mandate),
  time_end = 300, increment = 1.0
)

R

# create intervention_type column for plotting
output_masks$intervention_type <- "mask mandate"
output_baseline$intervention_type <- "baseline"
output <- rbind(output_masks, output_baseline)

ggplot(data = output[output$compartment == "infectious", ]) +
  aes(
    x = time,
    y = value,
    color = intervention_type,
    linetype = intervention_type
  ) +
  stat_summary(
    fun = "sum",
    geom = "line",
    linewidth = 1
  ) +
  scale_y_continuous(
    labels = scales::comma
  ) +
  labs(
    x = "Simulation time (days)",
    y = "Individuals"
  ) +
  theme_bw(
    base_size = 15
  ) +
  geom_vline(
    xintercept = c(mask_mandate$time_begin, mask_mandate$time_end),
    colour = "black",
    linetype = "dashed",
    linewidth = 0.2
  ) +
  annotate(
    geom = "text",
    label = "Mask mandate",
    colour = "black",
    x = (mask_mandate$time_end - mask_mandate$time_begin) / 2 +
      mask_mandate$time_begin,
    y = 10,
    angle = 0,
    vjust = "outward"
  )

Intervention types

There are two intervention types for model_default_cpp(). Rate interventions on model parameters (transmissibillity $\beta$, infectiousness_rate $\sigma$ and recovery_rate $\gamma$) and contact matrix reductions contacts.

To implement both contact and rate interventions in the same simulation they must be passed as a list e.g. intervention = list(transmissibility = mask_mandate, contacts = close_schools). But if there are multiple interventions that target contact rates, these must be passed as one contacts input. See the vignette on modelling overlapping interventions for more detail.

Pharmaceutical interventions

Pharmaceutical interventions (PIs) are measures such as vaccination and mass treatment programs. In the previous section, we assumed that interventions reduced the value of parameter values while the intervention was in place. In the case of vaccination, we assume that after the intervention individuals are no longer susceptible and should be classified into a different disease state. Therefore, we specify the rate at which individuals are vaccinated and track the number of vaccinated individuals through time.

The diagram below shows the SEIRV model implemented using model_default_cpp() where susceptible individuals are vaccinated and then move to the $V$ class.

The equations describing this model are as follows:

\[ \begin{aligned} \frac{dS_i}{dt} & = - \beta S_i \sum_j C_{i,j} I_j -\nu_{t} S_i \\ \frac{dE_i}{dt} &= \beta S_i\sum_j C_{i,j} I_j - \alpha E_i \\ \frac{dI_i}{dt} &= \alpha E_i - \gamma I_i \\ \frac{dR_i}{dt} &=\gamma I_i \\ \frac{dV_i}{dt} & =\nu_{i,t} S_i\\ \end{aligned} \] Individuals are vaccinated at an age group ($i$) specific time dependent ($t$) vaccination rate ($\nu_{i,t}$). The SEIR components of these equations are described in the tutorial simulating transmission.

To explore the effect of vaccination we need to create a vaccination object to pass as an input into model_default_cpp() that includes an age groups specific vaccination rate nu and age group specific start and end times of the vaccination program (time_begin and time_end).

Here we will assume all age groups are vaccinated at the same rate 0.01 and that the vaccination program starts on day 40 and is in place for 150 days.

R

# prepare a vaccination object
vaccinate <- vaccination(
  name = "vaccinate all",
  time_begin = matrix(40, nrow(contact_matrix)),
  time_end = matrix(40 + 150, nrow(contact_matrix)),
  nu = matrix(c(0.01, 0.01, 0.01))
)

We pass our vaccination object using vaccination = vaccinate:

R

output_vaccinate <- model_default_cpp(
  population = uk_population,
  transmissibility = 2.7 / 5.5,
  infectiousness_rate = 1.0 / 4.0,
  recovery_rate = 1.0 / 5.5,
  vaccination = vaccinate,
  time_end = 300, increment = 1.0
)

Compare interventions

Plot the three interventions vaccination, school closure and mask mandate and the baseline simulation on one plot. Which intervention reduces the peak number of infectious individuals the most?

Output

R

# create intervention_type column for plotting
output_vaccinate$intervention_type <- "vaccination"
output <- rbind(output_school, output_masks, output_vaccinate, output_baseline)

ggplot(data = output[output$compartment == "infectious", ]) +
  aes(
    x = time,
    y = value,
    color = intervention_type,
    linetype = intervention_type
  ) +
  stat_summary(
    fun = "sum",
    geom = "line",
    linewidth = 1
  ) +
  scale_y_continuous(
    labels = scales::comma
  ) +
  labs(
    x = "Simulation time (days)",
    y = "Individuals"
  ) +
  theme_bw(
    base_size = 15
  )

From the plot we see that the peak number of total number of infectious individuals when vaccination is in place is much lower compared to school closures and mask wearing interventions.

Summary

Different types of intervention can be implemented using mathematical modelling. Modelling interventions requires assumptions of which model parameters are affected (e.g. contact matrices, transmissibility), by what magnitude and and what times in the simulation of an outbreak.

The next step is to quantify the effect of an interventions. If you are interested in learning how to compare interventions, please complete the tutorial Comparing public health outcomes of interventions.

Key Points

The effect of NPIs can be modelled as reducing contact rates between age groups or reducing the transmissibility of infection
Vaccination can be modelled by assuming individuals move to a different disease state $V$

Content from Comparing public health outcomes of interventions

Last updated on 2024-09-24 | Edit this page

Estimated time: 75 minutes

Overview

Questions

How can I quantify the effect of an intervention?

Objectives

Compare intervention scenarios

Prerequisites

Complete tutorials Simulating transmission and Modelling interventions

Learners should familiarise themselves with following concept dependencies before working through this tutorial:

Outbreak response : Intervention types.

Introduction

In this tutorial we will compare intervention scenarios against each other. To quantify the effect of the intervention we need to compare our intervention scenario to a counter factual scenario. The counter factual is the scenario in which nothing changes, often referred to as the ‘do nothing’ scenario. The counter factual scenario may include no interventions, or if we are investigating the potential impact of an additional intervention in the later stages of an outbreak there may be existing interventions in place.

We must also decide what our outcome of interest is to make comparisons between intervention and counter factual scenarios. The outcome of interest can be:

a model outcome, e.g. number of infections or hospitalisations,
a metric such as the epidemic peak time or size,
a measure that uses the model outcomes such as QALY/DALYs.

Instructor Note

In this tutorial we introduce the concept of the counter factual and how to compare scenarios (counter factual versus intervention) against each other.

Vacamole model

The Vacamole model is a deterministic model based on a system of ODEs in Ainslie et al. 2022 to describe the effect of vaccination on COVID-19 dynamics. The model consists of 11 compartments, individuals are classed as one of the following:

susceptible, $S$,
partial vaccination ($V_1$), fully vaccination ($V_2$),
exposed, $E$ and exposed while vaccinated, $E_V$,
infectious, $I$ and infectious while vaccinated, $I_V$,
hospitalised, $H$ and hospitalised while vaccinated, $H_V$,
dead, $D$,
recovered, $R$.

The diagram below describes the flow of individuals through the different compartments.

Running a counterfactual scenario using the Vacamole model

Run the model with the default parameter values for the UK population assuming that :

1 in a million individual are infectious (and not vaccinated) at the start of the simulation
The contact matrix for the United Kingdom has age bins:
- age between 0 and 20 years,
- age between 20 and 40,
- 40 years and over.
There is no vaccination scheme in place

Using the output, plot the number of deaths through time

Vaccination code

To run the model with no vaccination in place we can either create two vaccination objects (one for each dose) using vaccination() with the time start, time end and vaccination rate all set to 0, or we can use the no_vaccination() function to create a vaccination object for two doses with all values set to 0.

R

no_vaccination <- no_vaccination(population = uk_population, doses = 2)

HINT : Running the model with default parameter values

We can run the Vacamole model with default parameter values by just specifying the population object and number of time steps to run the model for:

R

output <- model_vacamole_cpp(
  population = uk_population,
  vaccination = no_vaccination,
  time_end = 300
)

SOLUTION

Run the model

R

polymod <- socialmixr::polymod
contact_data <- socialmixr::contact_matrix(
  survey = polymod,
  countries = "United Kingdom",
  age.limits = c(0, 20, 40),
  symmetric = TRUE
)

OUTPUT

Using POLYMOD social contact data. To cite this in a publication, use the 'get_citation()' function

OUTPUT

Removing participants that have contacts without age information. To change this behaviour, set the 'missing.contact.age' option

R

# prepare contact matrix
contact_matrix <- t(contact_data$matrix)

# extract demography vector
demography_vector <- contact_data$demography$population
names(demography_vector) <- rownames(contact_matrix)

# prepare initial conditions
initial_i <- 1e-6

initial_conditions <- c(
  S = 1 - initial_i,
  V1 = 0, V2 = 0,
  E = 0, EV = 0,
  I = initial_i, IV = 0,
  H = 0, HV = 0, D = 0, R = 0
)

initial_conditions <- rbind(
  initial_conditions,
  initial_conditions,
  initial_conditions
)
rownames(initial_conditions) <- rownames(contact_matrix)

# prepare population object
uk_population <- population(
  name = "UK",
  contact_matrix = contact_matrix,
  demography_vector = demography_vector,
  initial_conditions = initial_conditions
)

ERROR

Error in population(name = "UK", contact_matrix = contact_matrix, demography_vector = demography_vector, : could not find function "population"

R

no_vaccination <- no_vaccination(population = uk_population, doses = 2)

ERROR

Error in no_vaccination(population = uk_population, doses = 2): could not find function "no_vaccination"

R

# run model
output <- model_vacamole_cpp(
  population = uk_population,
  vaccination = no_vaccination,
  time_end = 300
)

ERROR

Error in model_vacamole_cpp(population = uk_population, vaccination = no_vaccination, : could not find function "model_vacamole_cpp"

Plot the number of deaths through time

R

ggplot(output[output$compartment == "dead", ]) +
  geom_line(
    aes(time, value, colour = demography_group),
    linewidth = 1
  ) +
  scale_colour_brewer(
    palette = "Dark2",
    labels = rownames(contact_matrix),
    name = "Age group"
  ) +
  scale_y_continuous(
    labels = scales::comma
  ) +
  labs(
    x = "Simulation time (days)",
    y = "Individuals"
  ) +
  theme(
    legend.position = "top"
  ) +
  theme_bw(
    base_size = 15
  )

ERROR

Error in ggplot(output[output$compartment == "dead", ]): could not find function "ggplot"

Comparing scenarios

Coming soon

Challenge : Ebola outbreak analysis

Coming soon

Key Points

The counter factual scenario must be defined to make comparisons

Probability (\(p\))	Expected change
\(p < 0.05\)	Increasing
\(0.05 \leq p< 0.4\)	Likely increasing
\(0.4 \leq p< 0.6\)	Stable
\(0.6 \leq p < 0.95\)	Likely decreasing
\(0.95 \leq p \leq 1\)	Decreasing

Overview

Questions

Objectives

Prerequisites

Why to use R packages for Outbreak analytics?

Questions to think about

Example: Quantify transmission

R

First, get your case data

R

OUTPUT

Then, set the generation time

R

Let’s calculate the reproduction number!

R

OUTPUT

R

Is this \(Rt\) estimation biased?

The problem!

What can we do?

How can I start?

Key Points

Overview

Questions

Objectives

Prerequisites

Introduction

Definitions

R

R

R

Find a Generation time

From time periods to probability distributions.

Serial interval

Give me a hint

Show me the solution

Instructor Note

Extract epidemiological parameters

R

OUTPUT

OUTPUT

R

OUTPUT

OUTPUT

R

OUTPUT

OUTPUT

CASE-INSENSITIVE

R

OUTPUT

OUTPUT

Ebola’s incubation periods

Give me a hint

Show me the solution

R

Select a single distribution

R

OUTPUT

OUTPUT

R

OUTPUT

OUTPUT

How does single_epidist works?

R

OUTPUT

R

OUTPUT

Ebola’s incubation period

Give me a hint

R

Show me the solution

R

Extract the summary statistics

R

OUTPUT

How to get the sd and other nested elements?

Give me a hint

R

Show me the solution

R

How does `single_epidist` works?

How to get the `sd` and other nested elements?

Plug-in `{epiparameter}` to `{EpiNow2}`