---
title: "Age-stratified hospitalisation and death risks"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Age-stratified hospitalisation and death risks}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

::: {.alert .alert-info}
If you are unfamiliar with the {simulist} package or the `sim_linelist()` function [Get Started vignette](simulist.html) is a great place to start.
:::

This vignette will describe how to simulate a line list with age-stratified (or age-dependent) hospitalisation and death risks. 

There are three risks in {simulist}, specifically in the `sim_linelist()` and `sim_outbreak()` functions, that relate to hospitalisations and deaths:

1. Hospitalisation risk (`hosp_risk`)
2. Death risk in hospitals (`hosp_death_risk`)
3. Death risk outside of hospitals (`non_hosp_death_risk`)

The default for `sim_linelist()` is a `0.2` (or 20%) hospitalisation risk for individuals infected, `0.5` (or 50%) death risk for hospitalised individuals, and `0.05` (or 5%) death risk for infected individuals outside of hospitals. These default risks are applied to all age groups in the populations. 

However, it may be unrealistic to assume that the probability an infected person is admitted to hospital or dies is independent of their age. For many diseases, young children or elderly individuals are at higher risk of being hospitalised and having adverse outcomes.

The `sim_linelist()` function can accommodate age-stratified risks by accepting a `<data.frame>` instead of a single risk for the entire population.

```{r setup}
library(simulist)
library(epiparameter)
```

Here is an example that uses the default hospitalisation and death risks (inside and outside of hospital). First we load the requested delay distributions using the {epiparameter} package.

```{r, read-delay-dists}
contact_distribution <- epiparameter(
  disease = "COVID-19",
  epi_name = "contact distribution",
  prob_distribution = create_prob_distribution(
    prob_distribution = "pois",
    prob_distribution_params = c(mean = 2)
  )
)

infectious_period <- epiparameter(
  disease = "COVID-19",
  epi_name = "infectious period",
  prob_distribution = create_prob_distribution(
    prob_distribution = "gamma",
    prob_distribution_params = c(shape = 1, scale = 1)
  )
)

# get onset to hospital admission from {epiparameter} database
onset_to_hosp <- epiparameter_db(
  disease = "COVID-19",
  epi_name = "onset to hospitalisation",
  single_epiparameter = TRUE
)

# get onset to death from {epiparameter} database
onset_to_death <- epiparameter_db(
  disease = "COVID-19",
  epi_name = "onset to death",
  single_epiparameter = TRUE
)
```

We set the seed to ensure we have the same output each time the vignette is rendered. When using {simulist}, setting the seed is not required unless you need to simulate the same line list multiple times.

```{r, set-seed}
set.seed(1)
```

## Population-wide risks

Simulate a line list with population-wide default risks: 

* hospitalisation risk: `0.2`
* death risk in hospitals `0.5`
* death risk outside of hospitals `0.05`

```{r, sim-linelist}
linelist <- sim_linelist(
  contact_distribution = contact_distribution,
  infectious_period = infectious_period,
  prob_infection = 0.5,
  onset_to_hosp = onset_to_hosp,
  onset_to_death = onset_to_death
)

# first 6 rows of linelist
head(linelist)
```

We can run another simulation and change the hospitalisation and death risks, inside and outside the hospital, still applied to the entire population. In this scenario the probability of being hospitalised if infected is higher, but the mortality risk for both hospitalised and non-hospitalised groups is lower.

```{r, sim-linelist-diff-risks}
linelist <- sim_linelist(
  contact_distribution = contact_distribution,
  infectious_period = infectious_period,
  prob_infection = 0.5,
  onset_to_hosp = onset_to_hosp,
  onset_to_death = onset_to_death,
  hosp_risk = 0.4,
  hosp_death_risk = 0.2,
  non_hosp_death_risk = 0.01
)

head(linelist)
```

## Age-stratified hospitalisation and death risks

To define age-stratified risks, we must create a table (`<data.frame>`) which contains the lower limits of the age groups and their respective risks.

In this example the hospitalisation risk will be:

* `0.2` (or 20%) for those aged 80 years or older
* `0.1` (or 10%) for those younger than 5 years
* `0.05` (or 5%) for the remaining age group

The oldest age group stops at the upper age range given in the `population_age` argument. The default upper age range in 90, so in this example the oldest age bracket will be 80-90 (inclusive). The minimum age of each age group is inclusive, and the maximum age of each age group is exclusive, except the oldest age group which is inclusive of the minimum and maximum age. In this example the first age group is the first element of each vector, so the minimum age is 1, maximum age is four (as the next age group starts at five), and the hospitalisation risk for that group is 0.1. Each age group forms a row in the table.

```{r, make-hosp-risk-df}
age_dep_hosp_risk <- data.frame(
  age_limit = c(1, 5, 80),
  risk = c(0.1, 0.05, 0.2)
)
age_dep_hosp_risk
```

```{r, sim-age-strat-linelist}
linelist <- sim_linelist(
  contact_distribution = contact_distribution,
  infectious_period = infectious_period,
  prob_infection = 0.5,
  onset_to_hosp = onset_to_hosp,
  onset_to_death = onset_to_death,
  hosp_risk = age_dep_hosp_risk
)

head(linelist)
```

::: {.alert .alert-warning}

The minimum age of the youngest age group must match the age range specified in the `population_age` argument of `sim_linelist()`, and the largest age limit in the risk `<data.frame>` must not be older than the upper age range. If these conditions are not met the function will error.

If the age-stratified risk table does not match the default (`c(1, 90)`), the `population_age` argument will need to be set to match.

:::

For example, the default age range of the population is 1 to 90 (inclusive). In our example above, the lowest age group started at 1 and the oldest age group stopped at 90. This matches the default `population_age = c(1, 90)`. However, see here that if the lower age limit exceeds the age range the function will not run.

```{r, sim-age_strat_linelist-error, error=TRUE}
age_dep_hosp_risk <- data.frame(
  age_limit = c(1, 5, 95),
  risk = c(0.1, 0.05, 0.2)
)
age_dep_hosp_risk

linelist <- sim_linelist(
  contact_distribution = contact_distribution,
  infectious_period = infectious_period,
  prob_infection = 0.5,
  onset_to_hosp = onset_to_hosp,
  onset_to_death = onset_to_death,
  hosp_risk = age_dep_hosp_risk
)
```

In order to make this code run with the age-stratified hospitalisation risk given, the `population_age` can be adjusted. Here the oldest age bracket is now 95 to 100 (`[95, 100]`).

```{r, sim-age_strat-linelist-diff-age-range}
age_dep_hosp_risk <- data.frame(
  age_limit = c(1, 5, 95),
  risk = c(0.1, 0.05, 0.2)
)

linelist <- sim_linelist(
  contact_distribution = contact_distribution,
  infectious_period = infectious_period,
  prob_infection = 0.5,
  onset_to_hosp = onset_to_hosp,
  onset_to_death = onset_to_death,
  hosp_risk = age_dep_hosp_risk,
  population_age = c(1, 100)
)

head(linelist)
```

Exactly the same method of age-stratified risks applies to death risks. First create the `<data.frame>` with the age groups and their respective, in this case, death risks, and then supply this to either the `hosp_death_risk` or `non_hosp_death_risk` arguments, to define the death risks in and outside the hospital, respectively, or both.

Here are a couple of examples:

```{r, sim-age_strat_linelist-hosp-death-risk}
age_dep_hosp_death_risk <- data.frame(
  age_limit = c(1, 5, 80),
  risk = c(0.3, 0.1, 0.6)
)
age_dep_hosp_death_risk

linelist <- sim_linelist(
  contact_distribution = contact_distribution,
  infectious_period = infectious_period,
  prob_infection = 0.5,
  onset_to_hosp = onset_to_hosp,
  onset_to_death = onset_to_death,
  hosp_death_risk = age_dep_hosp_death_risk
)
```

```{r, sim-age_strat_linelist-non-hosp-death-risk}
age_dep_non_hosp_death_risk <- data.frame(
  age_limit = c(1, 5, 80),
  risk = c(0.1, 0.05, 0.1)
)
age_dep_non_hosp_death_risk

linelist <- sim_linelist(
  contact_distribution = contact_distribution,
  infectious_period = infectious_period,
  prob_infection = 0.5,
  onset_to_hosp = onset_to_hosp,
  onset_to_death = onset_to_death,
  non_hosp_death_risk = age_dep_non_hosp_death_risk
)
```

Up until now these age-stratified tables have only been supplied to one risk. However, these can be supplied in the same simulation. In this case the hospitalisation risk, and death risks inside and outside of hospital, are all specified.

```{r, sim-age-strat-linelist-all}
age_dep_hosp_risk <- data.frame(
  age_limit = c(1, 5, 80),
  risk = c(0.1, 0.05, 0.2)
)
age_dep_hosp_death_risk <- data.frame(
  age_limit = c(1, 5, 80),
  risk = c(0.3, 0.1, 0.6)
)
age_dep_non_hosp_death_risk <- data.frame(
  age_limit = c(1, 5, 80),
  risk = c(0.1, 0.05, 0.1)
)

linelist <- sim_linelist(
  contact_distribution = contact_distribution,
  infectious_period = infectious_period,
  prob_infection = 0.5,
  onset_to_hosp = onset_to_hosp,
  onset_to_death = onset_to_death,
  hosp_risk = age_dep_hosp_risk,
  hosp_death_risk = age_dep_hosp_death_risk,
  non_hosp_death_risk = age_dep_non_hosp_death_risk
)

head(linelist)
```