epidemics_odin.Rmd

---
title: "Getting started with epidemic scenario modelling components"
output:
  bookdown::html_vignette2:
    fig_caption: yes
    code_folding: show
pkgdown:
  as_is: true
bibliography: references.json
link-citations: true
vignette: >
  %\VignetteIndexEntry{Getting started with epidemic scenario modelling components}
  %\VignetteEncoding{UTF-8}
  %\VignetteEngine{knitr::rmarkdown}
editor_options: 
  chunk_output_type: console
---

This initial vignette shows to get started with using the _epidemics_ package.

Further vignettes include guidance on ["Modelling the implementation of vaccination regimes"](modelling_vaccination.html), as well as on["Modelling non-pharmaceutical interventions (NPIs) to reduce social contacts"](modelling_interventions.html) and ["Modelling multiple overlapping NPIs"](modelling_multiple_interventions.html).

There is also guidance available on specific models in the model library, such as the [Vacamole model developed by RIVM, the Dutch Institute for Public Health](model_vacamole.html).

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

```{r setup}
# library(epidemics)
devtools::load_all()
library(dplyr)
library(ggplot2)
```

## Prepare population and initial conditions

Prepare population and contact data.

::: {.alert .alert-info}
### Note on social contacts data {-}

**Note that** the social contacts matrices provided by the [_socialmixr_](https://CRAN.R-project.org/package=socialmixr) package follow a format wherein the matrix $M_{ij}$ represents [contacts from group $i$ to group $j$](https://epiforecasts.io/socialmixr/articles/socialmixr.html#usage).

However, epidemic models traditionally adopt the notation that $M_{ij}$ defines contacts to $i$ from $j$ [@wallinga2006].

$q M_{ij} / n_i$ then defines the probability of infection, where $q$ is a scaling factor dependent on $R_0$ (or another measure of infection transmissibility), and $n_i$ is the population proportion of group $i$.
The ODEs in _epidemics_ also follow this convention.

For consistency with this notation, social contact matrices from _socialmixr_ need to be transposed (using `t()`) before they are used with _epidemics_.
:::

```{r}
# load contact and population data from socialmixr::polymod
polymod <- socialmixr::polymod
contact_data <- socialmixr::contact_matrix(
  polymod,
  countries = "United Kingdom",
  age.limits = c(0, 20, 40),
  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)

# view contact matrix and demography
contact_matrix

demography_vector
```

Prepare initial conditions for each age group.

```{r}
# initial conditions
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 <- rbind(
  initial_conditions,
  initial_conditions,
  initial_conditions
)

# assign rownames for clarity
rownames(initial_conditions) <- rownames(contact_matrix)

# view initial conditions
initial_conditions
```

Prepare a population as a `population` class object.

```{r}
uk_population <- population(
  name = "UK",
  contact_matrix = contact_matrix,
  demography_vector = demography_vector,
  initial_conditions = initial_conditions
)

uk_population
```

## Run epidemic model

```{r}
# run an epidemic model using `epidemic`
output <- model_default(
  population = uk_population,
  time_end = 600, increment = 1.0
)

output_odin <- model_default_odin(
  population = uk_population,
  time_end = 600, increment = 1.0
)

library(data.table)
all.equal(output, output_odin, tolerance = 0.3) # FALSE
all.equal(output, output_odin, tolerance = 0.4) # TRUE
```

## Prepare data and visualise infections

Plot epidemic over time, showing only the number of individuals in the exposed and infected compartments.

```{r class.source = 'fold-hide'}
# plot figure of epidemic curve
filter(output, compartment %in% c("exposed", "infectious")) %>%
  ggplot(
    aes(
      x = time,
      y = value,
      col = demography_group,
      linetype = compartment
    )
  ) +
  geom_line() +
  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() +
  theme(
    legend.position = "top"
  ) +
  labs(
    x = "Simulation time (days)",
    linetype = "Compartment",
    y = "Individuals"
  )
```

## References