Skip to content

add compilation of edits #126

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
avallecam merged 3 commits into from
Nov 22, 2025
Merged

add compilation of edits #126

Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
92d48e8
add compilation of edits
avallecam Nov 22, 2025
dccb05d
fix lintr
avallecam Nov 22, 2025
3996d0c
remove intervention outputs + add namespace
avallecam Nov 22, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
144 changes: 80 additions & 64 deletions episodes/modelling-interventions.Rmd
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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}
Expand Down Expand Up @@ -53,64 +52,99 @@ 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 <- epidemics::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:

```{r intervention}
close_schools <- intervention(
close_schools <- epidemics::intervention(
name = "School closure",
type = "contacts",
time_begin = 50,
Expand All @@ -126,41 +160,28 @@ 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(
output_school <- epidemics::model_default(
# population
population = uk_population,
# rate
Expand All @@ -175,24 +196,13 @@ 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 <- base::rbind(output_school, output_baseline)

output %>%
filter(compartment == "infectious") %>%
Expand Down Expand Up @@ -224,11 +234,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.

Expand All @@ -237,7 +248,7 @@ We expect that mask wearing will reduce an individual's infectiousness, based on
We create an intervention object with `type = rate` and `reduction = 0.161`. Using parameters adapted from [Li et al. 2020](https://doi.org/10.1371/journal.pone.0237691) we have proportion wearing masks = coverage $\times$ availability = $0.54 \times 0.525 = 0.2835$ and proportion reduction in transmission rate = $0.575$. Therefore, $\theta = 0.2835 \times 0.575 = 0.163$. We assume that the mask wearing mandate starts at day 40 and continue to be in place for 200 days.

```{r masks}
mask_mandate <- intervention(
mask_mandate <- epidemics::intervention(
name = "mask mandate",
type = "rate",
time_begin = 40,
Expand All @@ -249,7 +260,7 @@ mask_mandate <- intervention(
To implement this intervention on the transmission rate $\beta$, we specify `intervention = list(transmission_rate = mask_mandate)`.

```{r output_masks}
output_masks <- model_default(
output_masks <- epidemics::model_default(
# population
population = uk_population,
# rate
Expand All @@ -268,7 +279,7 @@ output_masks <- model_default(
# create intervention_type column for plotting
output_masks$intervention_type <- "mask mandate"
output_baseline$intervention_type <- "baseline"
output <- rbind(output_masks, output_baseline)
output <- base::rbind(output_masks, output_baseline)

output %>%
filter(compartment == "infectious") %>%
Expand Down Expand Up @@ -370,18 +381,18 @@ Here we will assume all age groups are vaccinated at the same rate 0.01 and that

```{r vaccinate}
# prepare a vaccination object
vaccinate <- vaccination(
vaccinate <- epidemics::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))
)
```

We pass our vaccination object into the model using the argument `vaccination = vaccinate`:

```{r output_vaccinate}
output_vaccinate <- model_default(
output_vaccinate <- epidemics::model_default(
# population
population = uk_population,
# rate
Expand All @@ -408,7 +419,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 <- base::rbind(
output_school,
output_masks,
output_vaccinate,
output_baseline
)

output %>%
filter(compartment == "infectious") %>%
Expand Down
Loading