add compilation of edits

- create baseline model section
- change object names when using socialmixr
- expand the initial conditions vector code
- create matrix with simpler code
- print internvetion class objects
- reduce one level to subtitles of interventions
- expand binding of data frame when comparing all internvetions

Author: avallecam

episodes/modelling-interventions.Rmd

@@ -2,7 +2,6 @@
 title: 'Modelling interventions'
 teaching: 45 # teaching time in minutes
 exercises: 30 # exercise time in minutes
-
 ---
 
 ```{r setup, echo= FALSE, message = FALSE, warning = FALSE}
@@ -53,58 +52,93 @@ In this tutorial different types of intervention and how they can be modelled ar
 
 ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
 
-## Non-pharmaceutical interventions
+## A baseline model
 
-[Non-pharmaceutical interventions](../learners/reference.md#NPIs) (NPIs) are measures put in place to reduce transmission that do not include the administration of drugs or vaccinations. NPIs aim at reducing contacts between infectious and susceptible individuals by closure of schools and workplaces, and other measures to prevent the spread of the disease,  for example, washing hands and wearing masks.
+We will investigate the effect of interventions on a COVID-19 outbreak using an SEIR model (`model_default()` in the R package `{epidemics}`). To be able to see the effect of our intervention, we will run a baseline variant of the model, i.e,  without intervention.
 
-We will investigate the effect of interventions on a COVID-19 outbreak using an SEIR model (`model_default()` in the R package `{epidemics}`). The SEIR model divides the population into four compartments: Susceptible (S), Exposed (E), Infectious (I), and Recovered (R). We will set the following parameters for our model: $R_0 = 2.7$ (basic reproduction number), latent period or pre-infectious period $= 4$ days, and the infectious period $= 5.5$ days (parameters adapted from [Davies et al. (2020)](https://doi.org/10.1016/S2468-2667(20)30133-X)). We adopt a contact matrix with age bins 0-18, 18-65, 65 years and older using `{socialmixr}`, and assume that one in every 1 million individuals in each age group is infectious at the start of the epidemic.
+The SEIR model divides the population into four compartments: Susceptible (S), Exposed (E), Infectious (I), and Recovered (R). We will set the following parameters for our model: $R_0 = 2.7$ (basic reproduction number), latent period or pre-infectious period $= 4$ days, and the infectious period $= 5.5$ days (parameters adapted from [Davies et al. (2020)](https://doi.org/10.1016/S2468-2667(20)30133-X)). We adopt a contact matrix with age bins 0-18, 18-65, 65 years and older using `{socialmixr}`, and assume that one in every 1 million individuals in each age group is infectious at the start of the epidemic.
 
 ```{r model_setup, echo = TRUE, message = FALSE}
-polymod <- socialmixr::polymod
-contact_data <- socialmixr::contact_matrix(
-  polymod,
+# load survey data
+survey_data <- socialmixr::polymod
+
+# generate contact matrix
+cm_results <- socialmixr::contact_matrix(
+  survey = survey_data,
   countries = "United Kingdom",
   age.limits = c(0, 15, 65),
   symmetric = TRUE
 )
 
-# prepare contact matrix
-contact_matrix <- t(contact_data$matrix)
+# transpose contact matrix
+cm_matrix <- t(cm_results$matrix)
 
 # prepare the demography vector
-demography_vector <- contact_data$demography$population
-names(demography_vector) <- rownames(contact_matrix)
+demography_vector <- cm_results$demography$population
+names(demography_vector) <- rownames(cm_matrix)
 
 # initial conditions: one in every 1 million is infected
 initial_i <- 1e-6
 initial_conditions <- c(
-  S = 1 - initial_i, E = 0, I = initial_i, R = 0, V = 0
+  S = 1 - initial_i, 
+  E = 0, 
+  I = initial_i, 
+  R = 0, 
+  V = 0
 )
 
 # build for all age groups
-initial_conditions <- matrix(
-  rep(initial_conditions, dim(contact_matrix)[1]),
-  ncol = 5, byrow = TRUE
+initial_conditions <- base::rbind(
+  initial_conditions,
+  initial_conditions,
+  initial_conditions
 )
-rownames(initial_conditions) <- rownames(contact_matrix)
+rownames(initial_conditions) <- rownames(cm_matrix)
 
 # prepare the population to model as affected by the epidemic
 uk_population <- epidemics::population(
   name = "UK",
-  contact_matrix = contact_matrix,
+  contact_matrix = cm_matrix,
   demography_vector = demography_vector,
   initial_conditions = initial_conditions
 )
 ```
 
-#### Effect of school closures on COVID-19 spread
+We run the model with a transmission rate $= 2.7/5.5$ (remember that [transmission rate = $R_0$* recovery rate](../episodes/simulating-transmission.md#the-basic-reproduction-number-r_0)), infectiousness rate $1/= 4$ and the recovery rate $= 1/5.5$ as follows:
+
+```{r, echo = TRUE, message = FALSE}
+# time periods
+preinfectious_period <- 4.0
+infectious_period <- 5.5
+basic_reproduction <- 2.7
+
+# rates
+infectiousness_rate <- 1.0 / preinfectious_period
+recovery_rate <- 1.0 / infectious_period
+transmission_rate <- basic_reproduction * recovery_rate
+
+# run baseline simulation with no intervention
+output_baseline <- model_default(
+  population = uk_population,
+  transmission_rate = transmission_rate,
+  infectiousness_rate = infectiousness_rate,
+  recovery_rate = recovery_rate,
+  time_end = 300, increment = 1.0
+)
+```
+
+## Non-pharmaceutical interventions
+
+[Non-pharmaceutical interventions](../learners/reference.md#NPIs) (NPIs) are measures put in place to reduce transmission that do not include the administration of drugs or vaccinations. NPIs aim at reducing contacts between infectious and susceptible individuals by closure of schools and workplaces, and other measures to prevent the spread of the disease,  for example, washing hands and wearing masks.
+
+### Effect of school closures on COVID-19 spread
 
 The first NPI we will consider is the effect of school closures on reducing the number of individuals infected with COVID-19 over time. We assume that a school closure will reduce the frequency of contacts within and between different age groups. Based on empirical studies, we assume that school closures will reduce the contacts between school-aged children (aged 0-15) by 50%, and will cause a small reduction (1%) in the contacts between adults (aged 15 and over). 
 
 To include an intervention in our model we must create an `intervention` object. The inputs are the name of the intervention (`name`), the type of intervention (`contacts` or `rate`), the start time (`time_begin`), the end time (`time_end`) and the reduction (`reduction`). The values of the reduction matrix are specified in the same order as the age groups in the contact matrix. 
 
 ```{r}
-rownames(contact_matrix)
+rownames(cm_matrix)
 ```
 
 Therefore, we specify ` reduction = matrix(c(0.5, 0.01, 0.01))`. We assume that the school closures start on day 50 and continue to be in place for a further 100 days. Therefore our intervention object is: 
@@ -117,6 +151,8 @@ close_schools <- intervention(
   time_end = 50 + 100,
   reduction = matrix(c(0.5, 0.01, 0.01))
 )
+
+close_schools
 ```
 
 ::::::::::::::::::::::::::::::::::::: callout
@@ -126,38 +162,25 @@ In `{epidemics}`, the contact matrix is scaled down by proportions for the perio
 
 ```{r echo = FALSE}
 reduction <- matrix(c(0.5, 0.1))
-contact_matrix_example <- matrix(c(1, 1, 1, 1), nrow = 2)
-contact_matrix_example
+cm_matrix_example <- matrix(c(1, 1, 1, 1), nrow = 2)
+cm_matrix_example
 ```
 
 If the reduction is 50% in group 1 and 10% in group 2, the contact matrix during the intervention will be:
 
 ```{r echo = FALSE}
-contact_matrix_example[1, ] <- contact_matrix_example[1, ] * (1 - reduction[1])
-contact_matrix_example[, 1] <- contact_matrix_example[, 1] * (1 - reduction[1])
-contact_matrix_example[2, ] <- contact_matrix_example[2, ] * (1 - reduction[2])
-contact_matrix_example[, 2] <- contact_matrix_example[, 2] * (1 - reduction[2])
-contact_matrix_example
+cm_matrix_example[1, ] <- cm_matrix_example[1, ] * (1 - reduction[1])
+cm_matrix_example[, 1] <- cm_matrix_example[, 1] * (1 - reduction[1])
+cm_matrix_example[2, ] <- cm_matrix_example[2, ] * (1 - reduction[2])
+cm_matrix_example[, 2] <- cm_matrix_example[, 2] * (1 - reduction[2])
+cm_matrix_example
 ```
 
 The contacts within group 1 are reduced by 50% twice to accommodate for a 50% reduction in outgoing and incoming contacts ($1\times 0.5 \times 0.5 = 0.25$). Similarly, the contacts within group 2 are reduced by 10% twice. The contacts between group 1 and group 2 are reduced by 50% and then by 10% ($1 \times 0.5 \times 0.9= 0.45$). 
 
 ::::::::::::::::::::::::::::::::::::::::::::::::
 
-Using transmission rate $= 2.7/5.5$ (remember that [transmission rate = $R_0$/ infectious period](../episodes/simulating-transmission.md#the-basic-reproduction-number-r_0)), infectiousness rate $1/= 4$ and the recovery rate $= 1/5.5$ we run the model with ` intervention = list(contacts = close_schools)` as follows:
-
-```{r}
-# time periods
-preinfectious_period <- 4.0
-infectious_period <- 5.5
-basic_reproduction <- 2.7
-
-# rates
-infectiousness_rate <- 1.0 / preinfectious_period
-recovery_rate <- 1.0 / infectious_period
-transmission_rate <- basic_reproduction / infectious_period
-```
-
+We run the model with ` intervention = list(contacts = close_schools)` as follows:
 
 ```{r school}
 output_school <- model_default(
@@ -175,25 +198,14 @@ output_school <- model_default(
 ```
 
 
-To be able to see the effect of our intervention, we also run a baseline variant of the model, i.e,  without intervention, combine the two outputs into one data frame, and then plot the output. Here we plot the total number of infectious individuals in all age groups using `ggplot2::stat_summary()` function:
+To observe the effect of our intervention, we will combine the baseline and intervention outputs into a single data frame and then plot the results. Here we plot the total number of infectious individuals in all age groups using `ggplot2::stat_summary()` function:
 
 ```{r baseline, echo = TRUE, fig.width = 10}
-# run baseline simulation with no intervention
-output_baseline <- model_default(
-  population = uk_population,
-  transmission_rate = transmission_rate,
-  infectiousness_rate = infectiousness_rate,
-  recovery_rate = recovery_rate,
-  time_end = 300, increment = 1.0
-)
-
 # create intervention_type column for plotting
 output_school$intervention_type <- "school closure"
 output_baseline$intervention_type <- "baseline"
 output <- rbind(output_school, output_baseline)
 
-library(tidyverse)
-
 output %>%
   filter(compartment == "infectious") %>%
   ggplot() +
@@ -224,11 +236,12 @@ output %>%
     y = "Individuals"
   )
 ```
+
 We can see that with the intervention in place, the infection still spreads through the population and hence accumulation of immunity contributes to the eventual peak-and-decline. However, the peak number of infectious individuals is smaller (green dashed line) than the baseline with no intervention in place (red solid line), showing a reduction in the absolute number of cases.
 
 
 
-#### Effect of mask wearing on COVID-19 spread
+### Effect of mask wearing on COVID-19 spread
 
 We can also model the effect of other NPIs by reducing the value of the relevant parameters. For example, investigating the effect of mask wearing on the number of individuals infected with COVID-19 over time. 
 
@@ -244,6 +257,8 @@ mask_mandate <- intervention(
   time_end = 40 + 200,
   reduction = 0.163
 )
+
+mask_mandate
 ```
 
 To implement this intervention on the transmission rate $\beta$, we specify `intervention = list(transmission_rate = mask_mandate)`.
@@ -372,10 +387,12 @@ Here we will assume all age groups are vaccinated at the same rate 0.01 and that
 # prepare a vaccination object
 vaccinate <- vaccination(
   name = "vaccinate all",
-  time_begin = matrix(40, nrow(contact_matrix)),
-  time_end = matrix(40 + 150, nrow(contact_matrix)),
+  time_begin = matrix(40, nrow(cm_matrix)),
+  time_end = matrix(40 + 150, nrow(cm_matrix)),
   nu = matrix(c(0.01, 0.01, 0.01))
 )
+
+vaccinate
 ```
 
 We pass our vaccination object into the model using the argument `vaccination = vaccinate`:
@@ -408,7 +425,12 @@ Plot the three interventions vaccination, school closure and mask mandate and th
 ```{r plot_vaccinate, echo = TRUE, message = FALSE, fig.width = 10}
 # create intervention_type column for plotting
 output_vaccinate$intervention_type <- "vaccination"
-output <- rbind(output_school, output_masks, output_vaccinate, output_baseline)
+output <- rbind(
+  output_school,
+  output_masks,
+  output_vaccinate,
+  output_baseline
+)
 
 output %>%
   filter(compartment == "infectious") %>%