---
title: "Exploring uncertainty in pandemic costs due to parameter uncertainty"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Exploring uncertainty in pandemic costs due to parameter uncertainty}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

# set seed for local reproducibility
set.seed(1)
```

This vignette shows how to explore the effect of model parameter uncertainty on epidemic costs using _daedalus_.

Here we explore the effect of uncertainty in $R_0$, which is a common use case in the initial stages of a novel pandemic when estimates of $R_0$ are less robust.

We shall model uncertainty in a scenario of a SARS-2004 like infection outbreak in a country chosen at random, Norway.

We focus on infection parameter uncertainty, but a similar approach can be applied to country characteristics.

```{r setup}
# load {daedalus} and helper packages
library(daedalus)

library(data.table)
library(ggplot2)
library(ggdist)
```

## General approach

The general approach is:

- Create multiple `<daedalus_infection>` objects with distinct values of parameters;

- Run `daedalus()` on each distinct `<daedalus_infection>` object;

- Get costs from each model output object.


## Drawing parameters from a distribution

First, we access the SARS 2004 $R_0$ provided in _daedalus_, and using this as a mean we draw 100 samples from a normal distribution with a standard deviation of 0.1.

```{r get_sars_r0}
# get SARS-1 R0
sars_2004 <- daedalus_infection("sars_cov_1")
sars_2004_r0 <- get_data(sars_2004, "r0")

# draw samples from a normal distribution
r0_samples <- rnorm(100, sars_2004_r0, sd = 0.1)
```

## Create infection objects

Next we create `<daedalus_infection>` objects to hold the $R_0$ values.

**Note that** future versions of _daedalus_ will include better support for parameter uncertainty, such as that provided by the [_distributional_ package](https://cran.r-project.org/package=distributional).

```{r make_inf_objs}
# make a list of infection objects
infection_list <- lapply(
  r0_samples, function(x) daedalus_infection("sars_cov_1", r0 = x)
)

# View an infection object
infection_list[[1]]
```

## Run DAEDALUS model for each infection object

We run the model for Norway with each infection object, representing different values of $R_0$.

We assume that there is no vaccine investment in advance of the outbreak. Vaccination with a pathogen specific vaccine is thus assumed to begin 1 year after the epidemic start date, similar to the situation for Covid-19.

We also assume for this example that this outbreak is unmitigated and allowed to become an epidemic.

The model runs for 600 days or a little over 1.5 years.

**Note that** running `daedalus()` iteratively may take some time, as the function checks the inputs and internally prepares parameters each time, as well as preparing the output data. Future package versions aim to streamline these steps and provide a stripped down version of the function optimised for scenario modelling for parameter fitting or with parameter uncertainty.

Users can run the model iteratively over the $R_0$ samples in parallel using frameworks such as [the _furrr_ package](https://cran.r-project.org/package=furrr).

```{r uncertainty_run_model}
# run daedalus()
output_list <- lapply(
  infection_list, daedalus,
  country = "Norway"
)
```

## Get epidemic costs

We can get the total epidemic costs --- a combination of life years lost, educational losses leading to lost future earnings, and economic losses due to worker illness --- using the helper function `get_costs()` and passing the option `summarise_as = "total"`.

```{r}
# get costs
costs <- vapply(
  output_list, get_costs, numeric(1), "total"
)
```

We can plot the total costs to view the distribution using the [_ggdist_ package](https://cran.r-project.org/package=ggdist) to visualise the distribution.

```{r plot_total_costs_dist}
# Use {ggplot2} to visualise the output
# `stat_histinterval()` is from {ggdist}
ggplot() +
  stat_histinterval(
    aes(costs),
    fill = "steelblue"
  ) +
  scale_x_continuous(
    labels = scales::comma_format()
  ) +
  theme_tidybayes() +
  labs(
    x = "Overall epidemic cost (million $)",
    y = "Density"
  )
```