---
title: "Getting started with epidist"
description: "A quick start guide to using the epidist R package"
output:
bookdown::html_vignette2:
fig_caption: yes
code_folding: show
number_sections: true
pkgdown:
as_is: true
# csl: https://raw.githubusercontent.com/citation-style-language/styles/master/apa-numeric-superscript-brackets.csl
link-citations: true
vignette: >
%\VignetteIndexEntry{Getting started with epidist}
%\VignetteEngine{knitr::rmarkdown}
%\VignetteEncoding{UTF-8}
bibliography: references.bib
---
`epidist` is a toolkit for flexibly estimating epidemiological delays built on top of `brms` (Bayesian Regression Models using Stan).
`brms` provides a powerful framework for Bayesian modelling with an accessible R interface to Stan.
By building on `brms`, `epidist` inherits the ability to work within the broader `brms` ecosystem, allowing users to leverage existing tools for model diagnostics, posterior predictive checks, and model comparison while addressing the specific challenges of delay distribution estimation.
See the `vignette("faq")` for more details on the tools available in the `brms` ecosystem.
In this vignette, we will give a quick start guide to using the `epidist` package.
To get started we will introduce some of the key concepts in delay distribution estimation, and then simulate some data delay data from a stochastic outbreak that includes common biases.
Using this simulated data we will then show how to use the `epidist` package to estimate a distribution using a simple model and a model that accounts for some of the common issues in delay distribution estimation.
We will then compare the output of these models to the true delay distribution used to simulate the data again using `epidst` tools.
# Key concepts in delay distribution estimation
In epidemiology, we often need to understand the time between key events - what we call "delays".
Think about things like:
* incubation period (how long between getting infected and showing symptoms),
* serial interval (time between when one person shows symptoms and when someone they infected shows symptoms), and
* generation interval (time between when one person gets infected and when they infect someone else).
We can think of all these as the time between a "primary event" and a "secondary event".
The tricky bit? Getting accurate estimates of these delays from real-world data is a challenge. The two main challenges we typically face are:
1. interval censoring (we often only know events happened within a time window, not the exact time), and
2. right truncation (we might miss observing later events if our observation period ends).
Don't worry if these terms sound a bit technical! In Section \@ref(data), we'll walk through what these issues look like by simulating the kind of data you might see during an outbreak.
Then in Section \@ref(fit), we'll show how `epidist` helps you estimate delay distributions accurately by accounting for these issues.
::: {.alert .alert-primary}
For those interested in the technical details, `epidist` implements models following best practices in the field.
Check out @park2024estimating for a methodological overview and @charniga2024best for a practical checklist designed for applied users.
We also recommend the [nowcasting and forecasting infectious disease dynamics](https://nfidd.github.io/nfidd/) course for more hands on learning.
:::
# Setup {#setup}
To run this vignette yourself, as well as the `epidist` package, you will need the following packages:
``` r
library(epidist)
library(ggplot2)
library(dplyr)
```
# Simulating data {#data}
We simulate data from an outbreak setting where the primary event is symptom onset and the secondary event is case notification.
We assume that both events are dates and so we do not know precise event times.
This is typically the most common setting for delay distribution estimation.
This is a simplified version of the more complete setting that `epidist` supports where events can have different censoring intervals.
We also assume that we are observing a sample of cases during the outbreak.
This means that our data is both interval censored for each event and truncated for the secondary event.
We first assume that the reporting delay is lognormal with a mean log of 1.6 and a log standard deviation of 0.5.
Here we use the `add_summaries()` function to add the mean and sd to the `data.frame`.
``` r
secondary_dist <- data.frame(mu = 1.6, sigma = 0.5) |>
add_summaries(family = "lognormal")
secondary_dist
#> mu sigma mean sd
#> 1 1.6 0.5 5.612521 2.991139
```
We can visualise the delay distribution these parameters describe with `plot()` (Figure \@ref(fig:lognormal)).
(ref:lognormal) The lognormal distribution is skewed to the right. Long delay times still have some probability.
``` r
plot(secondary_dist, type = "delay", max_delay = 30)
```
We then simulate the stochastic outbreak in continuous time, sample a reporting delay for each event and finally observe the events (i.e restrict them to be dates).
We assume that the outbreak has a growth rate of 0.2, that we observe the outbreak for 25 days, and that the first case was infected on the 1st of February 2024.
``` r
growth_rate <- 0.2
obs_time <- 25
outbreak_start_date <- as.Date("2024-02-01")
```
Click to expand for simulation details
First, we use the [Gillespie algorithm](https://en.wikipedia.org/wiki/Gillespie_algorithm) to generate infectious disease outbreak data from a stochastic compartmental model.
``` r
outbreak <- simulate_gillespie(r = growth_rate, seed = 101)
```
`outbreak` is a `data.frame` with the two columns: `case` and `ptime`.
Here `ptime` is a numeric column giving the time of infection.
In reality, it is more common to receive primary event times as a date rather than a numeric.
``` r
head(outbreak)
#> case ptime
#> 1 1 0.04884052
#> 2 2 0.06583120
#> 3 3 0.21857827
#> 4 4 0.24963421
#> 5 5 0.30133392
#> 6 6 0.31425010
```
To generate secondary events, we will use a lognormal distribution (Figure \@ref(fig:lognormal)) for the delay between primary and secondary events:
``` r
obs <- simulate_secondary(
outbreak,
dist = rlnorm,
meanlog = secondary_dist[["mu"]],
sdlog = secondary_dist[["sigma"]]
)
```
`obs` is now a `data.frame` with further columns for `delay` and `stime`.
The secondary event time is simply the primary event time plus the delay:
``` r
all(obs$ptime + obs$delay == obs$stime)
#> [1] TRUE
```
If we were to receive the complete data `obs` as above then it would be simple to accurately estimate the delay distribution.
However, in reality, during an outbreak we almost never receive the data as above.
First, the times of primary and secondary events will usually be censored.
This means that rather than exact event times, we observe event times within an interval.
Here we suppose that the interval is daily, meaning that only the date of the primary or secondary event, not the exact event time, is reported (Figure \@ref(fig:cens)):
``` r
obs_cens <- simulate_dates(
obs,
outbreak_start_date = outbreak_start_date,
obs_time = obs_time
)
head(obs_cens)
#> case pdate_lwr pdate_upr sdate_lwr sdate_upr obs_date
#> 1 1 2024-02-01 2024-02-02 2024-02-05 2024-02-06 2024-02-26
#> 2 2 2024-02-01 2024-02-02 2024-02-09 2024-02-10 2024-02-26
#> 3 3 2024-02-01 2024-02-02 2024-02-07 2024-02-08 2024-02-26
#> 4 4 2024-02-01 2024-02-02 2024-02-09 2024-02-10 2024-02-26
#> 5 5 2024-02-01 2024-02-02 2024-02-04 2024-02-05 2024-02-26
#> 6 6 2024-02-01 2024-02-02 2024-02-08 2024-02-09 2024-02-26
```
`simulate_dates()` floors each event time to its reporting window, daily by default, offsets it from the date the outbreak started, and records the date of the last observation.
As a final step we rename the columns to show that `epidist` does not require particular names.
``` r
obs_data <- transmute(
obs_cens,
id = case,
symptom_onset = pdate_lwr,
case_notification = sdate_lwr,
obs_date = obs_date
)
```
Issues not considered
Another issue, which `epidist` currently does not account for, is that sometimes only the secondary event might be observed, and not the primary event.
For example, symptom onset may be reported, but start of infection unknown.
Discarding events of this type leads to what are called ascertainment biases.
Whereas each case is equally likely to appear in the sample above, under ascertainment bias some cases are more likely to appear in the data than others.
`linelist_data` is the data we will fit to (Figure \@ref(fig:linelist)).
(ref:linelist) The primary and secondary event windows of each observed case, ordered by the date of symptom onset. The dashed line is the observation date, after which no secondary event is seen. These are the 200 sampled cases, drawn from across the whole outbreak rather than the first 500 cases of Figure \@ref(fig:cens).
``` r
plot_events(linelist_data, obs_time = max(linelist_data$obs_date))
```
Click to expand for code to create the observed data histogram
``` r
delay_histogram <- plot_delays(
list(
"Censored retrospective data" = linelist_all,
"Censored, truncated,\nsampled data" = linelist_data
),
reference = c(
mu = secondary_dist[["mu"]],
sigma = secondary_dist[["sigma"]]
)
)
```
