---
title: "Estimating severity from data"
author: "Pablo N Perez-Guzman (materials adapted from Robert Verity)"
format: html
editor: visual
---

## Hierarchy of outcomes

As you saw from the slides, it is useful to imagine a *hierarchy of outcomes* when thinking about disease severity (**Figure 1**). At the base of this hierarchy are all (*true*) infections, some of which may lead to asymptomatic episodes that remain undetectable in clinical data. Above this, we have mild cases that are detectable via active surveillance (i.e., going out and measuring people), but are likely to be under-ascertained by passive surveillance (i.e., observing cases presenting at clinics) because they are mild enough that people tend to forgo healthcare. Above this, we have severe cases that warrant intensive medical intervention. In highly resourced settings we can think of these as cases that present to hospitals, but in resource-poor settings many of these severe cases may be missing from clinical data. Finally, at the top of the hierarchy, we have deaths. These tend to occur in clinical settings and tend to be the best recorded of all outcomes. However, in resource-poor settings there may be issues with vital registration, along with large numbers of community deaths that make hospital deaths a substantial underestimate of total deaths.

::: {style="text-align: center;"}
![<b>Figure 1:</b> The hierarchy of outcomes](images/pyramid.png){width="70%" height="70%"}
:::

We can quantify severity by comparing any one stage of this hierarchy against another. For example, the case fatality ratio (CFR) is the proportion of cases that result in death. Likewise, the infection fatality ratio (IFR) is the proportion of infections that result in death. We can navigate this hierarchy by multiplying values together, for example:

<br>

$$
\begin{array}{c}
\text{Infection} \\
\text{fatality ratio} \\
\text{(IFR)}
\end{array}
\hspace{5mm} = \hspace{5mm}
\begin{array}{c}
\text{Infection} \\
\text{hospitalisation ratio} \\
\text{(IHR)}
\end{array}
\hspace{5mm} \times \hspace{5mm}
\begin{array}{c}
\text{Hospitalisation} \\
\text{fatality ratio} \\
\text{(HFR)}
\end{array}
$$ Breaking the problem down into chunks like this allows us to use different sources of information to explore different parts of hierarchy before linking them together in an overall picture of the disease.

# Simple (crude) severity estimates

The simplest way to do a crude severity estimate $e(s)$ is to take the cumulative number of severe outcomes up to the present day $s$ and divide by cumulative cases. For example, to calculate the HFR we would do:

$$ HFR(s) = \frac{D(s)}{H(s)} $$

where $D(s)$ and $H(s)$ denote the cumulative deaths and hospitalisations up to time $s$.

::: {.callout-note icon="false"}
## Question 1: *What type of data would you need to have from a given outbreak, in order to make this calculation? What do you need to assume about those data, so that your HFR calculation is, indeed, an accurate severity estimate?*

There are several types of outbreak surveillance data that can be used to calculate severity. For this extremely simple calculation we do not anything other than the cumulative numbers up to the present day.

That said, two common outbreak surveillance data streams are aggregated time-series data and linelist data. Aggregated time-series data is basically a count of severity outcomes (hospitalisations and deaths, in our example) per unit-time. The reporting unit-time can be days, weeks or months (even years could make sense for chronic diseases), and the counts per unit-time can in turn correspond to incident (e.g., daily new cases) or cumulative outcomes up to that point in time.

Linelist data, on the other hand, is patient-level information on disease outcomes of interest. It will usually consist of a data set where each row corresponds to an individual patient, and dates (or number of days since the start of the outbreak) at which outcomes of interest occur, sometimes alongside other relevant clinical information.

This simple HFR definition assumes that 1) all deaths occur in hospital settings (i.e., all those infected who eventually die receive in-hospital care); and b) by this token, that all severe cases become hospitalisations that will either recover or die.
:::

Run the code in the next box to load our first example of outbreak data that we will use to calculate HFR.

```{r}
#| message: false
#| warning: false

# load necessary packages and support functions
source("support.R")
pacman::p_load(rio,            # import/export files
               tidyverse,      # data cleaning, wrangling, plotting, etc.
               drjacoby,       # flexible MCMC inference implementation
               patchwork,      # multi-panel plotting tools
               scales,         # labelling tools for plots
               matrixStats,    # stats on matrix object tools
               epitools,       # epi analysis tools
               cfr,            # tools to estimate disease severity
               kableExtra,     # tools to visualise data tables
               htmltools)      # tools for html-style visualisations

# read in raw data
data <- rio::import("incidence_data.csv")

# visualise data
data |>
  kbl() |>
  kable_styling(bootstrap_options = "striped", 
                full_width = F, position = "center") |>
  scroll_box(width = "900px", height = "200px")
```

::: {.callout-note icon="false"}
## Question 2: *What type of data is this and how would you then calculate the "true" HFR for this outbreak?*

This is a daily time-series of incident (new) infections, hospitalisations, recoveries and deaths. Assuming these indeed represent the *true* number of all outcomes observed per day of the outbreak (i.e., perfect surveillance system), we can divide the cumulative (total) number of deaths by cumulative hospitalisations by the end of the outbreak.
:::

Run the code below and check your intuition in the previous answer.

```{r}
#| message: false
#| warning: false
basic_hfr <- data %>% 
  select(day, deaths, hospitalisations) %>% 
  mutate(deaths = cumsum(deaths),
         hospitalisations = cumsum(hospitalisations)) %>% 
  mutate(HFR = deaths / hospitalisations) 

true_hfr <- tail(basic_hfr$HFR, 1)

ggplot(basic_hfr, aes(day, HFR)) +
  geom_point() +
  geom_hline(yintercept = true_hfr, col = "red", linetype = 2) +
  annotate("text", label = "True HFR", x = 25, y = 0.35) +
  labs(x = "Time", y = "HFR") +
  scale_y_continuous(limits = c(0, 0.5), expand = c(0, 0),
                     labels = scales::percent) +
  theme_minimal() +
  theme(axis.line = element_line())
```

::: {.callout-note icon="false"}
## Question 3: *Describe what you observe during this outbreak. What is happening with the HFR estimation during the early, middle and late stages of the outbreak?*

The outbreak has a total duration of 150 days. In the early and mid stages, there are several active (open) infections and hospitalisations for which we do not yet know the outcome. Our simple HFR calculation at these points is severely underestimating the *true* HFR, which we can only know until the outbreak is finished or at least until the very late stages, when enough cumulative hospitalisations and deaths have accrued.
:::

As we discussed during the slides, in a real-time outbreak situation data is often incomplete and there are reporting delays. Let's pretend we are now in the early stages of this outbreak.

In this scenario, the outbreak was detected 20 days after the first true infections occurred and, since, we have received data every 5 days and we are on day 30.

A simple way for gauging uncertainty in binomially-distributed observations, such as our simple HFR estimate, would be to calculate [binomial confidence intervals](https://en.wikipedia.org/wiki/Binomial_proportion_confidence_interval). Run the code below and compare what our simple estimator is telling us at each of the three data survey points we have so far (i.e., days 20, 25 and 30 of the outbreak); we will use a function built-in the package `epitools` to help with the calculation.

```{r}
#| message: false
#| warning: false
data <- rio::import("incidence_data_day30.csv")

n_hosp <- sum(data$hospitalisations)
n_death <- sum(data$deaths)
n_recov <- sum(data$recoveries)
n_open <- n_hosp - n_recov - n_death


basic_hfr <- data %>% 
  select(day, deaths, hospitalisations) %>% 
  mutate(deaths = cumsum(deaths),
         hospitalisations = cumsum(hospitalisations))

survey_times <- seq(20, 30, 5)

hfr_conf <- basic_hfr %>% 
  select(day, deaths, hospitalisations) %>% 
  filter(day %in% survey_times) %>% 
  mutate(epitools::binom.exact(x = deaths, n = hospitalisations)) %>% 
  mutate(mean = round(proportion, 3),
         lb = round(lower, 3),
         ub = round(upper, 3)) %>% 
  select(day, mean, lb, ub)


ggplot(hfr_conf) +
  geom_pointrange(aes(x = day, y = mean, ymin = lb, ymax = ub)) +
  geom_hline(yintercept = true_hfr, col = "red", linetype = 2) +
  annotate("text", label = "True HFR", x = 25, y = 0.35) +
  labs(x = "Time", y = "HFR") +
  scale_y_continuous(limits = c(0, 0.5), expand = c(0, 0),
                     labels = scales::percent) +
  theme_minimal() +
  theme(axis.line = element_line())
```

::: {.callout-note icon="false"}
## Question 4: *Describe what you observe.*

Once again, we see that the simple HFR calculation underestimates the *true* HFR, which for the sake of the example we are plotting here from the final outbreak data. This cannot be known at the start of an outbreak, however, and whilst the calculation of a binomial confidence interval adds some information on the data available to make our simple calculation, it does not do a good job at capturing the true HFR.
:::

Let's then look at another correction we can apply when all possible outcomes of a disease severity pathway are know. In our example of an outbreak of some disease causing acute illness, hospitalised patients can either recover or die. So, another way of calculating $HFR(s)$ could be

$$ HFR(s) = \frac{D(s)}{D(s) + R(s)}, $$

where $R(s)$ represents the cumulative number of hospitalised patients who have recovered up to time point $s$.

::: {.callout-note icon="false"}
## Question 5: *How do you think this new calculation will vary compared to the previous one where we divided D(s) by H(s))?*

The problems with the previous estimator were caused by cases being counted before they had transitioned into deaths. We can avoid this by only counting "closed" cases, i.e., those for which we know the final outcome. This is only possible if we have a measure of the total number of recoveries in the population, which may not always be true.
:::

Run the code below and check if your intuition was correct.

```{r}
#| message: false
#| warning: false

closed_hfr <- data %>% 
  select(day, deaths, hospitalisations, recoveries) %>% 
  filter(day %in% survey_times) %>% 
  mutate(deaths = cumsum(deaths),
         hospitalisations = cumsum(hospitalisations),
         recoveries = cumsum(recoveries)) %>% 
  mutate(epitools::binom.exact(x = deaths, 
                               n = recoveries + deaths)) %>% 
  mutate(mean = round(proportion, 3),
         lb = round(lower, 3),
         ub = round(upper, 3)) %>% 
  select(day, mean, lb, ub)


compare <- closed_hfr %>% 
  mutate(type = "Deaths / (Deaths + Recoveries)") %>% 
  rbind(., 
        hfr_conf %>% 
          mutate(type = "Deaths / Hospitalisations"))


ggplot(compare) +
  geom_pointrange(aes(x = day, y = mean, ymin = lb, ymax = ub, col = type)) +
  geom_hline(yintercept = true_hfr, col = "red", linetype = 2) +
  labs(x = "Time", y = "HFR") +
  scale_y_continuous(limits = c(0, 1), expand = c(0, 0),
                     labels = scales::percent) +
  theme_minimal() +
  theme(axis.line = element_line(),
        legend.position = c(0.7, 0.85),
        legend.title = element_blank())
```

::: {.callout-note icon="false"}
## Question 6: *What do you observe?*

The new estimate appears to be less biased, but has much greater uncertainty. This is because we are only making use of information in hospitalisations that have completed their journey to death or recovery, which in the early stages of an outbreak will be a small minority of cases.

There are also additional limitations of this approach that are worth noting. For one, if the average time to death is very different from the average time to recovery then estimates will be biased. It may also be difficult to obtain unbiased estimates of the actual number of recoveries as surveillance tends to be biased towards better recording of the more severe outcomes.
:::

# Accounting for delays

The simple estimates we have calculated so far demonstrate how we can make use of even minimal information, like aggregate numbers, to gauge severity. In some settings, that might be all we have. However, as we have seen, these estimates could be biased and unreliable, particularly during the early stages of an epidemic when most deaths have yet to occur.

Other more powerful methods that make use of open cases and explicitly deal with delays to outcomes are needed to avoid such biases. For this, let's look at the other type of data we mentioned, linelist data.

Run the code below to see some simulated linelist from our outbreak:

```{r}
linelist <- rio::import("simulated_linelist_day30.csv")
head(arrange(linelist, id), 20)


```

As you can see, each row represents an individual with a unique identifying number, their outcome and the delay (days) between hospitalisation and recovery or death. However, for most of the 508 hospitalised cases so far, the outcome is not yet known!

Run the code below to find out more about the number of hospitalised cases for which we have information on their outcomes and time to recovery or death.

```{r}
#| message: false
#| warning: false

n_recov <- sum(linelist$outcome == "recovery")
n_death <- sum(linelist$outcome == "death")
n_open <- nrow(linelist) - n_recov - n_death

delay_recov <- mean(linelist$delay_recovery, na.rm = TRUE)
delay_death <- mean(linelist$delay_death, na.rm = TRUE)


(ggplot(linelist) +
    geom_histogram(aes(delay_recovery)) +
    geom_vline(xintercept = delay_recov, linetype = 2, col = "red") +
    theme_minimal() +
    labs(x = "Delay to recovery", y = "n",
         title = paste("N =", n_recov, ", out of", nrow(linelist)))) +
  
  (ggplot(linelist) +
     geom_histogram(aes(delay_death)) +
     geom_vline(xintercept = delay_death, linetype = 2, col = "red") +
     theme_minimal() +
     labs(x = "Delay to death", y = "n",
          title = paste("N =", n_death, ", out of", nrow(linelist))
          ))

```

::: {.callout-note icon="false"}
## Question 7: What do you observe?

These distributions are quite noisy, especially deaths given that there have been relatively few "closed" hospitalisations to date. However, even with this limited data we can get a picture of the overall mean and spread. We also get a picture of the proportion of closed hospitalised cases that have resulted in death versus recovery.

The time from hospitalisation to recovery is quite variable, ranging from 0 to 28 days, with a mean around 12 days post-hospitalisation. Although there are fewer deaths, the observed time from hospitalisation to death ranged from 1 to 29 days, with a mean of 16 days post-hospitalisations.
:::

In our simulated data up to day 30 of the outbreak, 391 hospitalisations remain open. Estimates relying solely on closed cases may be biased. By dealing with censoring we can make use of the additional information from open cases as well - potentially leading to more accurate estimates of key parameters like the HFR. Two methods for doing this include the use of **parametric mixture models**, which we will come back to at the end of this practical, and the **Kaplan-Meier method** which we will discuss next.

Kaplan-Meier methods provide a **non-parametric way to account for censoring**. The full approach is described [here](https://doi.org/10.1093/aje/kwi230), but in short, we can graphically represent the number of cases that result in death or recovery using a Kaplan-Meier curve:

```{r}
#| message: false
#| warning: false
# get a vector of delays for closed recoveries and deaths
closed_recovery <- linelist %>% 
  filter(!is.na(delay_recovery)) %>% 
  pull(delay_recovery)

closed_death <- linelist %>% 
  filter(!is.na(delay_death)) %>% 
  pull(delay_death)

n <- nrow(linelist)


# calculate KM curve to recovery or death
df_km <- data.frame(day = 0:350) %>%
  mutate(
    n_recovered = 
      mapply(function(x) sum(closed_recovery <= x), day),
    n_dead = 
      mapply(function(x) sum(closed_death <= x), day),
    prop1 = 1 - n_recovered / n,
    prop2 = n_dead / n
  )

# plot KM curve
df_km %>% 
  select(day, prop1, prop2) %>% 
  pivot_longer(cols = -day) %>% 
  ggplot() + 
  geom_step(aes(x = day, y = value, col = name)) +
  labs(x = "Time from admission to outcome", y = "") +
  scale_y_continuous(limits = c(0, 1), 
                     breaks = seq(0, 1, 0.1),
                     expand = c(0, 0),
                     labels = scales::percent) +
  scale_x_continuous(expand = c(0, 0), limits = c(0, 30)) +
  scale_color_manual(values = c("chartreuse3", "purple"),
                     name = "", 
                     labels = c("Survivorship to\nrecovery", 
                                "\nOne minus\nsurvivorship to\ndeath")) +
  theme_minimal() +
  theme(axis.line = element_line(),
        legend.position = c(0.85, 0.5))
```

::: {.callout-note icon="false"}
## Question 8: What does the Y-axis and, more importantly, the gap between the two lines represent? If we had more data from later in the outbreak, do you think the lines would eventually intersect? What would such intersection represent?

The vertical gap between the lines indicates the proportion of cases without a recorded outcome. Once the epidemic is over and all cases have reached an outcome, the two lines will intersect. The point of intersection will be an unbiased estimate of the HFR. While the epidemic is ongoing, our estimate of the HFR plotted in the Y-axis must lie somewhere between the two lines.

Note that with this non-parametric method we can extrapolate and arrive at a point estimate of the HFR by making an implicit assumptioin that the relative probabilities of death versus recovery among open cases are the same as those observed in closed cases up to the time of the analysis. There is currently no known R package for performing this analysis, so we will not show an illustrative example here, but the approach is implemented in the [STATA package casefat](https://econpapers.repec.org/scripts/search.pf?ft=casefat).
:::

Line lists are extremely valuable when dealing with an emerging threat. They provide rich information at the individual level and can be used to account for right-censoring, leading to unbiased estimates of fatality rates and other proportions. However, we tend only to have line lists from passive surveillance, and so they tend to concentrate on the more severe end of the clinical spectrum.

Let's see what we can learn from more aggregate population-level time series data and how we can explicitly account for delays and censoring. We will once again use our daily incidence data, however we will now use the entire time series.

One important factor for this analysis is that we must have uninterrupted time series data, i.e. we cannot have any dates with missing data. If we do have missing data, then we should first impute those values before carrying out this analysis.

Run the code below to plot the entire time series of hospitalisations and deaths:

```{r}
#| message: false
#| warning: false

ts_data <- rio::import("incidence_data.csv") %>% 
  select(day, hospitalisations, deaths)

plot_data(ts_data, all_data = TRUE)

```

Looking at the incidence data, we can see a clear lag between the hospitalisations and deaths curve. The `cfr` package in R takes this lag into account and we can use it to estimate the HFR parametrically. That is, it **estimates the number of cases that are still open** based on the incidence data and an **assumed delay distribution between cases and deaths**. It then estimates the severity outcome of interest (HFR in our example) based on the estimated number of closed cases, which avoids some of the biases we saw earlier in our simple (crude) calculations.

Let's start with a delay distribution representing the probability that a hospitalised case dies at each time interval after symptom onset. In our example, we will assume a **parametric** form for this distribution - a [gamma distribution](https://en.wikipedia.org/wiki/Gamma_distribution). When doing this in a real-time outbreak response, we may not always have access to linelist data to estimate the **mean** and **standard deviation** of our gamma distribution, so we may have to make assumptions (e.g., based on knowledge from previous/similar outbreaks).

Here, we will get these parameters from our linelist data.

Run the code below to estimate a mean and standard deviation for our parametric model of the probability of death over time, and plot the resulting distribution:

```{r}
#| message: false
#| warning: false

discrete_gamma <- function(x, mean, sd) {
  pgamma(x, shape = mean^2 / sd^2, rate = mean / sd^2) -
    pgamma(x - 1, shape = mean^2 / sd^2, rate = mean / sd^2)
}


mean_delay <- mean(linelist$delay_death, na.rm = TRUE)

sd_delay <- sd(linelist$delay_death, na.rm = TRUE)


data.frame(x = 1:50) |>
  mutate(y = discrete_gamma(x, mean = mean_delay, sd = sd_delay)) |>
  ggplot() + theme_bw() +
  geom_bar(aes(x = x, y = y), stat = "identity") +
  xlab("Time (days)") + ylab("Probability of death")

```

We can now feed this delay distribution and our daily incidence data into the `cfr_static()` function, which would give us the following result:

```{r}
#| message: false
#| warning: false

# reformat data for cfr package
inc_df <- ts_data %>% 
  select(date = day,
         cases = hospitalisations,
         deaths) %>% 
  mutate(date = as.Date("2024-12-31") + date)

# estimate CFR accounting for delays
cfr_static(data = inc_df,
           delay_density = function(x) 
             discrete_gamma(x, mean = mean_delay, sd = sd_delay)) %>% 
  cbind(true_hfr)

```

In this example, we used our incomplete knowledge of the delay distribution to death from the linelist data we had up to day 30 of the outbreak. With these dat,a the `cfr_static` function estimated an HFR of 33.3% (95%CI 32.7-33.9), compared to our *true* estimate of 33.3% (i.e., based on total cumulative hospitalisations and deaths at the end of the outbreak). So, not bad at all!

Admittedly, the estimates are not different because the time-series we are employing now includes information past the peak of the outbreak, by which point there are few open hospitalised cases left. Let's see, however, how would this parametric method would perfomr if we compute HFR **throughout** the course of the outbreak with data available up to present time.

Run the code below to do this and compare the parametric estimate to our simple method based on cumulative counts:

```{r}
#| message: false
#| warning: false

day_start <- 20
# estimate CFR using cfr_static() from day 10 onward
df_cfr <- mapply(function(i) {
  cfr_static(data = inc_df[1:i,],
             delay_density = function(x) discrete_gamma(x, 
                                                        mean = mean_delay, 
                                                        sd = sd_delay)) %>% 
    suppressWarnings() %>% 
    suppressMessages()}, 
  day_start:nrow(inc_df), SIMPLIFY = FALSE) %>% 
  bind_rows() %>% 
  mutate(date = inc_df$date[day_start:nrow(inc_df)], .before = 1)

# crude estimate
hfr_df <- inc_df %>% 
  mutate(hfr_naive = cumsum(deaths) / cumsum(cases)) %>% 
  select(date, hfr_naive)

df_cfr %>% 
  left_join(hfr_df) %>% 
  ggplot(aes(x = date)) +
  geom_line(aes(y = hfr_naive, col = "Crude estimate")) +
  geom_line(aes(y = severity_estimate, col = "cfr_static()")) +
  geom_hline(yintercept = true_hfr, col = "red", linetype = 2) +
  geom_ribbon(aes(ymin = severity_low, ymax = severity_high), 
              fill = "green4", alpha = 0.2) +
  scale_color_manual(labels = c("cfr_static()", "Crude estimate"),
                     values = c("green4", "blue4")) +
  scale_y_continuous(limits = c(0, 1), expand = c(0, 0),
                     labels = scales::percent) +
  labs(x = "Date", y = "HFR") +
  theme_minimal() +
  theme(axis.line = element_line(),
        legend.title = element_blank(),
        legend.position = c(0.85, 0.85))
```

We can see that the `cfr_static()` function begins producing reliable estimates at an earlier point in the outbreak, compared to the simple estimate based on total counts. This is because the parametric model accounts for open cases and their expected delay distribution to a severe outcome.

Note that this parametric approach assumes we know the correct delay distribution, i.e., that the values of mean and SD are correct for this dataset. In our example, we know this because we have simulated the data! **In reality we will not know this distribution perfectly**, and will have to estimate it from data. In general, we are better off accounting for delays even if this distribution is imperfect rather than ignoring delays completely.

## Estimating delays and HFR together

Rather than assuming a known delay distribution, we can **estimate this distribution directly from the data**. If the time series of deaths looks like a lagged version of the time series of cases, then this lag can tell us the average time from case to death. Similarly, if the time series of deaths looks like a more spread out version of the time series of cases, then this tells us something about the variation in time from case to death.

We can combine these ideas into a mathematical model. That is, we can assume a simple parametric form for the delay distribution defined in terms of a mean and standard deviation. Although there are no out-of-the-box solutions for writing and fitting this mathematical model, there are many tools available to perform such type of Bayesian inference.

Here, we will use [drjacoby](https://mrc-ide.github.io/drjacoby/), an in-house package developed at Imperial to performs **Bayesian inference via Markov chain Mote Carlo (MCMC)**. The Bayesian approach requires us to define priors on the unknown parameters, which are distributions that capture what we believe about the parameters *before* we have seen any data. Alongside incorporating prior knowledge, the Bayesian framework allows us to:

-   **Estimate the HFR**, along with the delay distribution's **mean** and **standard deviation**
-   **Quantify uncertainty** around all parameter estimates

For this example, we will use the following weakly informative priors:

-   **HFR**: Uniform prior that is equally likely to take any value between 0 and 1
-   **Mean delay**: could take any positive value but is most likely to be between 0 and 90 days
-   **Standard deviation of the delay**: can take any positive value but is most likely to be between 0 and 30 days.

```{r}
#| message: false
#| warning: false

# function for convolving a vector of cases with a delay distribution
case_convolve <- function(x, delay) {
  ret <- rev(convolve(x = rev(x), y = delay, type = "open"))[1:length(x)]
  ret[ret < 0] <- 0
  return(ret)
}


# define parameters data.frame; we will use our estimated mean and SD delays as initial values
df_params <- define_params(name = "HFR", min = 0, max = 1, init = 0.1,
                           name = "delay_death_mean", min = 0, max = Inf, init = mean_delay,
                           name = "delay_death_sd", min = 0, max = Inf, init = sd_delay)

# define log-likelihood function
loglike <- function(params, data, misc) {
  
  # extract parameter values
  HFR <- params["HFR"]
  delay_death_mean <- params["delay_death_mean"]
  delay_death_sd <- params["delay_death_sd"]
  
  # get incidence of death by convolution
  n <- nrow(data)
  inc_death <- 
    HFR * case_convolve(x = data$cases,
                        delay = discrete_gamma(0:100, 
                                               mean = delay_death_mean, 
                                               sd = delay_death_sd))
  
  # calculate likelihood
  ret <- sum(dpois(data$deaths, lambda = inc_death[1:n], log = TRUE))
  
  # return
  return(ret)
}

# define log-prior function
logprior <- function(params, misc) {
  
  # extract parameter values
  HFR <- params["HFR"]
  delay_death_mean <- params["delay_death_mean"]
  delay_death_sd <- params["delay_death_sd"]
  
  # calculate log-prior
  ret <- dunif(HFR, min = 0, max = 1, log = TRUE) +
    dexp(delay_death_mean, rate = 1 / 16, log = TRUE) +
    dexp(delay_death_sd, rate = 1 / 8.5, log = TRUE)
  
  # return
  return(ret)
}

# run MCMC
set.seed(1)
mcmc <- run_mcmc(data = inc_df,
                 df_params = df_params,
                 loglike = loglike,
                 logprior = logprior,
                 burnin = 1e2,
                 samples = 1e3,
                 chains = 5,
                 silent = TRUE)


# get 95% CrIs
mcmc_draws <- mcmc$output |>
  filter(phase == "sampling") |>
  select(-chain, -phase, -iteration, -logprior, -loglikelihood)

mcmc_CrI <- apply(mcmc_draws, 2, function(x) quantile(x, probs = c(0.025, 0.5, 0.975))) |>
  t() |>
  as.data.frame() |>
  rename(lower_CrI = '2.5%',
         upper_CrI = '97.5%',
         median = '50%')

# show table of data
mcmc_CrI |>
  round(digits = 3) |>
  kbl() |>
  kable_styling(bootstrap_options = "striped", full_width = FALSE, position = "center")
```

There are many steps involved in running an MCMC and checking that results are reliable. We will not dwell on these steps, but briefly (in MCMC lingo) this involves assessing some convergence statistics and visual inspection of traceplots. Let's forget about these for now and instead focus on the actual results.

We obtained a mean and 95% credible interval estimates about the HFR, and also on the mean and standard deviation of the delay distribution to death after hospitalisation. Unlike the previous analysis, any uncertainty in the parameters of the delay distribution are taken into account in our estimate of the HFR. Results are similar to the previous analysis, but we have had to make fewer assumptions!

### Summary

-   Crude estimates can be useful, especially in situations where we have very little data to estimate severity from. Simple corrections to these can be easily made, by accounting for closed cases in our denominator.
-   Kaplan-Meier methods are useful to explicitly account for delay distributions to the outcomes of interest. Their application requires having access to linelist data, which may not always be possible.
-   The `cfr` package accounts for delays between cases and deaths and does not require linelist data. This parametric approach can provide more accurate severity estimates, by adjusting for open cases. It requires making assumptions about the delay distribution to outcome and likely parameters.
-   Bayesian inference, such as the MCMC we implemented with `drjacoby`, offers a flexible framework for jointly estimating the severity outcome and the delay distribution. It requires a bit more technical knowledge to implement as there are no out-of-the-box methods for it, but also offers great advantages such as incorporating uncertainty in both the severity outcome probability and delay distribution, and having to make less assumptions about the available data.
