---
title: "Modelling interventions that change infection parameters"
output:
  bookdown::html_vignette2:
    fig_caption: yes
    code_folding: show
pkgdown:
  as_is: true
vignette: >
  %\VignetteIndexEntry{Modelling interventions that change infection parameters}
  %\VignetteEncoding{UTF-8}
  %\VignetteEngine{knitr::rmarkdown}
editor_options: 
  chunk_output_type: console
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  message = FALSE,
  warning = FALSE,
  fig.width = 5,
  fig.height = 4,
  dpi = 300
)
```

::: {.alert .alert-warning}
**New to _epidemics_, or to modelling interventions?** It may help to read the ["Get started"](epidemics.html) first! See the ["Modelling a non-pharmaceutical intervention"](modelling_interventions.html) vignettes for a guide to modelling interventions on social contacts instead.
:::

```{r setup}
library(epidemics)
library(dplyr)
library(ggplot2)
```

## Prepare population and initial conditions

We prepare population and contact data from the U.K., with epidemiological compartments matching the default epidemic model (SEIR-V).

We assume that one in every million people has been infected and is infectious, translating to about 67 total infections for a U.K. population of 67 million.

The code for these steps is similar to that in the ["Getting started vignette"](epidemics.html) and is hidden here, although it can be expanded for reference.

```{r class.source = 'fold-hide'}
# 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, 65),
  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)
```

```{r class.source = 'fold-hide'}
# initial conditions
initial_i <- 1e-4
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)
```

```{r class.source = 'fold-hide'}
uk_population <- population(
  name = "UK",
  contact_matrix = contact_matrix,
  demography_vector = demography_vector,
  initial_conditions = initial_conditions
)
```

## Modelling an intervention on the transmission rate

We model an intervention on the transmission rate $\beta$ that reduces it by 10%.
This could represent interventions such as requiring people to wear masks that reduce transmission.

```{r}
# prepare an intervention that models mask mandates for ~3 months (100 days)
mask_mandate <- intervention(
  name = "mask mandate",
  type = "rate",
  time_begin = 60,
  time_end = 60 + 100,
  reduction = 0.1
)

# examine the intervention object
mask_mandate

# check the object
is_intervention(mask_mandate)

is_contacts_intervention(mask_mandate)

is_rate_intervention(mask_mandate)
```

We first run a baseline scenario --- no interventions are implemented to slow the spread of the epidemic --- and visualise the outcomes in terms of daily new infections.
We simulate an epidemic using `model_default()`, calling the default model outlined in the ["Get started vignette"](epidemics.html).

To examine the effect of a mask mandate, we simulate the epidemic for 200 days as we expect the intervention to spread disease incidence out over a longer period.

```{r}
# no intervention baseline scenario
data <- model_default(
  population = uk_population,
  time_end = 200, increment = 1.0
)

# with a mask mandate
data_masks <- model_default(
  population = uk_population,
  intervention = list(transmission_rate = mask_mandate),
  time_end = 200, increment = 1.0
)
```

```{r}
# get new infections in each scenario
data <- new_infections(data, by_group = TRUE)
data_masks <- new_infections(data_masks, by_group = TRUE)

# assign a scenario name to each scenario
data$scenario <- "baseline"
data_masks$scenario <- "masks"

# bind data together
data_combined <- bind_rows(data, data_masks)
```

We plot the data to examine the effect that implementing a mask mandate has on the daily number of new infections.

```{r class.source = 'fold-hide'}
ggplot(data_combined) +
  geom_line(
    aes(time, new_infections, col = demography_group, linetype = scenario)
  ) +
  coord_cartesian(
    expand = FALSE
  ) +
  annotate(
    geom = "rect",
    xmin = mask_mandate[["time_begin"]],
    xmax = mask_mandate[["time_end"]],
    ymin = 0, ymax = 150e3,
    fill = alpha("red", alpha = 0.2),
    lty = "dashed"
  ) +
  scale_y_continuous(
    labels = scales::comma
  ) +
  scale_linetype_manual(
    name = "Scenario",
    values = c(
      baseline = "dashed",
      masks = "solid"
    )
  ) +
  scale_colour_brewer(
    palette = "Dark2",
    name = "Age group"
  ) +
  expand_limits(
    y = c(0, 100e3)
  ) +
  coord_cartesian(
    expand = FALSE
  ) +
  theme_bw() +
  theme(
    legend.position = "top"
  ) +
  labs(
    x = "Simulation time (days)",
    linetype = "Compartment",
    y = "New infections"
  )
```