add episode on contact matrices

config.yaml

@@ -59,11 +59,13 @@ contact: '[email protected]'
 
 # Order of episodes in your lesson
 episodes:
+- contact-matrices.Rmd
 - simulating-transmission.Rmd
 - model-choices.Rmd
 - modelling-interventions.Rmd
 - compare-interventions.Rmd
 
+
 # Information for Learners
 learners:
 
@@ -78,6 +80,6 @@ profiles:
 # This space below is where custom yaml items (e.g. pinning
 # sandpaper and varnish versions) should live
 
+
 varnish: epiverse-trace/varnish@epiversetheme
-# this is carpentries/sandpaper#533 in our fork so we can keep it up to date with main
 sandpaper: epiverse-trace/sandpaper@patch-renv-github-bug

episodes/contact-matrices.Rmd

@@ -0,0 +1,271 @@
+---
+title: 'Contact matrices'
+teaching: 40
+exercises: 10
+---
+
+:::::::::::::::::::::::::::::::::::::: questions 
+
+- What is a contact matrix?
+- How are contact matrices estimated?
+- How are contact matrices used in epidemiological analysis?
+
+::::::::::::::::::::::::::::::::::::::::::::::::
+
+::::::::::::::::::::::::::::::::::::: objectives
+
+- Use the R package `socialmixr` to estimate a contact matrix
+- Understand the different types of analysis contact matrices can be used for 
+::::::::::::::::::::::::::::::::::::::::::::::::
+
+<!-- ::::::::::::::::::::::::::::::::::::: prereq -->
+
+<!-- ## Prerequisites -->
+
+
+
+<!-- ::::::::::::::::::::::::::::::::: -->
+
+
+## Introduction
+
+Some groups of individuals have more contacts than others; the average schoolchild has many more daily contact than the average elderly person, for example. This heterogeneity of contact patterns between different groups can affect disease transmission, because certain groups are more likely to transmit to others within that group, as well as to other groups. The rate at which individuals within and between groups make contact with others can be summarised in a contact matrix. In this tutorial we are going to learn how contact matrices can be used in different analyses and how the `{socialmixr}` package can be used to estimate contact matrices. 
+
+
+```{r,message=FALSE,warning=FALSE}
+library(socialmixr)
+```
+
+
+:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: instructor
+
+
+
+::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+
+## The contact matrix
+
+The basic contact matrix represents the amount of contact or mixing within and between different subgroups of a population. The subgroups are often age categories but can also be geographic areas or high/low risk groups. For example, a hypothetical contact matrix representing the average number of contacts per day between children and adults could be:
+
+$$
+\begin{bmatrix}
+2 & 2\\
+1 & 3 
+\end{bmatrix}
+$$
+In this example, we would use this to represent that children meet, on average, 2 other children and 2 adult per day (first row), and adults meet, on average, 1 child and 3 other adults per day (second row). We can use this kind of information to account for the role heterogeneity in contact plays in infectious disease transmission.
+
+
+## Using `socialmixr`
+
+Contact matrices are commonly estimated from studies that use diaries to record interactions. For example, the POLYMOD survey measured contact patterns in 8 European countries using data on the location and duration of contacts reported by the study participants [(Mossong et al. 2008)](https://doi.org/10.1371/journal.pmed.0050074).
+
+The R package `{socialmixr}` contains functions which can estimate contact matrices from POLYMOD and other surveys. We can load the POLYMOD survey data:
+
+
+```{r polymod_, echo = TRUE}
+polymod <- socialmixr::polymod
+```
+
+Then we can obtain the contact matrix for the age categories we want by specifying `age.limits`. 
+
+```{r polymod_uk, echo = TRUE}
+contact_data <- contact_matrix(
+  survey = polymod,
+  countries = "United Kingdom",
+  age.limits = c(0, 20, 40),
+  symmetric = TRUE
+)
+contact_data
+```
+
+
+
+**Note: although the contact matrix `contact_data$matrix` is not itself mathematically symmetric, it satisfies the condition that the total number of contacts of one group with another is the same as the reverse. In other words:
+`contact_data$matrix[j,i]*contact_data$demography$proportion[j] = contact_data$matrix[i,j]*contact_data$demography$proportion[i]`.
+For the mathematical explanation see [the corresponding section in the socialmixr documentation](https://epiforecasts.io/socialmixr/articles/socialmixr.html#symmetric-contact-matrices).**
+
+
+::::::::::::::::::::::::::::::::::::: callout
+### Why would a contact matrix be non-symmetric?
+
+One of the arguments we gave the function `contact_matrix()` is `symmetric=TRUE`. This ensures that the total number of contacts of age group 1 with age group 2 is the same as the total number of contacts of age group 2 and age group 1 (see the `socialmixr` [vignette](https://cran.r-project.org/web/packages/socialmixr/vignettes/socialmixr.html) for more detail). However, when contact matrices are estimated from surveys or other sources, the *reported* number of contacts may differ by age group resulting in a non-symmetric contact matrix because of biases in recall or reporting across different groups and/or uncertainty from using a limited sample of participants [(Prem et al 2021)](https://doi.org/10.1371/journal.pcbi.1009098). If `symmetric` is set to TRUE, the `contact_matrix()` function will internally use an average of reported contacts to ensure resulting total number of contacts are symmetric.
+
+::::::::::::::::::::::::::::::::::::::::::::::::
+
+The example above uses the POLYMOD survey. There are a number of surveys available in `socialmixr`, to list the available surveys use `list_surveys()`. To download a survey, we can use `get_survey()`
+
+```{r, message = FALSE, warning = FALSE}
+zambia_sa_survey <- get_survey("https://doi.org/10.5281/zenodo.3874675")
+```
+
+
+
+::::::::::::::::::::::::::::::::::::: challenge 
+
+## Zambia contact matrix
+
+After downloading the survey, generate a symmetric contact matrix for Zambia using the following age bins:
+
++ [0,15)
++ 20+
+
+:::::::::::::::::::::::: solution 
+
+```{r polymod_poland}
+contact_data_zambia <- contact_matrix(
+  survey = zambia_sa_survey,
+  age.limits = c(0, 20),
+  symmetric = TRUE
+)
+contact_data_zambia
+```
+:::::::::::::::::::::::::::::::::
+::::::::::::::::::::::::::::::::::::::::::::::::
+
+
+::::::::::::::::::::::::::::::::::::: callout
+## Synthetic contact matrices
+
+Contact matrices can be estimated from data obtained from diary (such as POLYMOD), survey or contact data, or synthetic ones can be used. [Prem et al. 2021](https://doi.org/10.1371/journal.pcbi.1009098) used the POLYMOD data within a Bayesian hierarchical model to project contact matrices for 177 other countries.
+
+::::::::::::::::::::::::::::::::::::::::::::::::
+
+
+
+
+## Analyses with contact matrices
+
+Contact matrices can be used in a wide range of epidemiological analyses, they can be used:
+
++ to calculate the basic reproduction number while accounting for different rates of contacts between age groups [(Funk et al. 2019)](https://doi.org/10.1186/s12916-019-1413-7),
++ to calculate final size of an epidemic, as in the R package `{finalsize}`,
++ to assess the impact of interventions finding the relative change between pre and post intervention contact matrices to calculate the relative difference in $R_0$ [(Jarvis et al. 2020)](https://doi.org/10.1186/s12916-020-01597-8),
++ and in mathematical models of tranmission within a population, to account for group specific contact patterns.
+
+
+However, all of these applications require us to perform some additional calculations using the contact matrix. Specifically, there are two main calculations we often need to do:
+
+1. **Convert contact matrix into expected number of secondary cases**
+
+If contacts vary between groups, then the average number of secondary cases won't be equal simply to the average number of contacts multiplied by the probability of transmission-per-contact. This is because the average amount of transmission in each generation of infection isn't just a matter of whom a group came into contact with; it's about whom *their contacts* subsequently come into contact with. The function `r_eff` in the package `{finalsize}` can perform this conversion, taking a contact matrix, demography and proportion susceptible and converting it into an estimate of the average number of secondary cases generated by a typical infectious individual (i.e. the effective reproduction number).
+
+2. **Convert contact matrix into contact rates**
+
+Whereas a contact matrix gives the average number of contacts that one groups makes with another, epidemic dynamics in different groups depend on the rate at which one group infects another. We therefore need to scale the rate of interaction between different groups (i.e. the number of contacts per unit time) to get the rate of transmission. However, we need to be careful that we are defining transmission to and from each group correctly in any model. Specifically, the entry $(i,j)$ in a mathematical model contact matrix represents contacts by group $i$ made with group $j$. But if we want to know the rate at which a group $i$ are getting infected, then we want to calculate the proportion of individuals infected in each group $j$ ($I_j/N_j$), multiplied by the number of contacts made with group $i$ ($C[i,j]$), multiplied by the transmission risk per contact ($\beta$).
+
+
+
+### In mathematical models
+
+Consider the SIR model where individuals are categorized as either susceptible $S$, infected but not yet infectious $E$, infectious $I$ or recovered $R$. The schematic below shows the processes which describe the flow of individuals between the disease states $S$, $I$ and $R$ and the key parameters for each process.
+
+```{r diagram, echo = FALSE, message = FALSE}
+DiagrammeR::grViz("digraph {
+
+  # graph statement
+  #################
+  graph [layout = dot,
+         rankdir = LR,
+         overlap = true,
+         fontsize = 10]
+
+  # nodes
+  #######
+  node [shape = square,
+       fixedsize = true
+       width = 1.3]
+
+       S
+       I
+       R
+
+  # edges
+  #######
+  S -> I [label = ' infection \n(transmission rate &beta;)']
+  I -> R [label = ' recovery \n(recovery rate &gamma;)']
+
+}")
+```
+
+The [differential equations](../learners/reference.md#ordinary) below describe how individuals move from one state to another [(Bjørnstad et al. 2020)](https://doi.org/10.1038/s41592-020-0822-z).
+
+
+$$
+\begin{aligned}
+\frac{dS}{dt} & = - \beta S I /N \\
+\frac{dI}{dt} &= \beta S I /N - \gamma I \\
+\frac{dR}{dt} &=\gamma I \\
+\end{aligned}
+$$
+To add age structure to our model, we need to add additional equations for the infection states $S$, $I$ and $R$ for each age group $i$. If we want to assume that there is heterogeneity in contacts between age groups then we must adapt the transmission term $\beta SI$ to include the contact matrix $C$ as follows :
+
+$$ \beta S_i \sum_j C_{i,j} I_j/N_j. $$ 
+
+Susceptible individuals in age group $i$ become infected dependent on their rate of contact with individuals in each age group. For each disease state ($S$, $E$, $I$ and $R$) and age group ($i$), we have a differential equation describing the rate of change with respect to time.  
+
+$$
+\begin{aligned}
+\frac{dS_i}{dt} & = - \beta S_i \sum_j C_{i,j} I_j/N_j \\
+\frac{dI_i}{dt} &= \beta S_i\sum_j C_{i,j} I_j/N_j - \gamma I_i \\
+\frac{dR_i}{dt} &=\gamma I_i \\
+\end{aligned}
+$$
+
+::::::::::::::::::::::::::::::::::::: callout
+##  Normalising the contact matrix to ensure the correct value of $R_0$
+
+When simulating an epidemic, we often want to ensure that the average number of secondary cases generated by a typical infectious individual (i.e. $R_0$) is consistent with known values for the pathogen we're analysing. In the above model, we scale the contact matrix by the $\beta$ to convert the raw interaction data into a transmission rate. But how do we define the value of $\beta$ to ensure a certain value of $R_0$?
+
+Rather than just using the raw number of contacts, we can instead normalise the contact matrix to make it easier to work in terms of $R_0$. Normalisation means converting to a value to be between 0 and 1. In particular, we normalise the matrix by scaling it so that if we were to calculate the average number of secondary cases based on this normalised matrix, the result would be 1 (in mathematical terms, we are scaling the matrix so the largest eigenvalue is 1). This transformation scales the entries but preserves their relative values.
+
+In the case of the above model, we want to define $\beta  C_{i,j}$ so that the model has a specified valued of $R_0$. $C[i,j]$ is defined as contacts to $i$ from $j$, which is equivalent to `contact_data$matrix[j,i]` so the first step is to transpose the contact data matrix (`contact_data$matrix`) so the row/column entries are now in the form $C[i,j]$. Then we normalise the matrix $C$ so the maximum eigenvalue is one and call this matrix $C_normalised$. Because the rate of recovery is $\gamma$, individuals will be infectious on average for $1/\gamma$ days. So $\beta$ is calculated from the scaling factor and the value of $\gamma$  (i.e. mathematically we have the dominant eigenvalue of the matrix $R_0 \times C_{normalised}$ is $\beta / \gamma$).
+
+```{r}
+contact_matrix <- t(contact_data$matrix)
+scaling_factor <- 1 / max(eigen(contact_matrix)$values)
+normalised_matrix <- contact_matrix * scaling_factor
+```
+
+As a result, if we multiply the scaled matrix by $R_0$, then converting to the number of expected secondary cases would give us $R_0$, as required.
+
+
+```{r}
+infectious_period <- 7.0
+basic_reproduction <- 1.46
+transmission_rate <- basic_reproduction * scaling_factor / infectious_period
+# check the dominant eigenvalue of R0 x C_normalised is R0
+max(eigen(basic_reproduction * normalised_matrix)$values)
+```
+
+
+::::::::::::::::::::::::::::::::::::::::::::::::
+
+
+::::::::::::::::::::::::::::::::::::: callout
+### Check the dimension of $\beta$
+
+In the SIR model without age structure the rate of contact is part of the transmission rate $\beta$, where as in the age-structured model we have separated out the rate of contact, hence the transmission rate $\beta$ in the age structured model will have a different value.
+
+::::::::::::::::::::::::::::::::::::::::::::::::
+
+We can use contact matrices from `socialmixr` with mathematical models in the R package `{epidemics}`. See the tutorial [Simulating transmission](../episodes/simulating-transmission.md) for examples and an introduction to `epidemics`.
+
+
+### Contact groups
+
+In the example above the dimension of the contact matrix will be the same as the number of age groups i.e. if there are 3 age groups then the contact matrix will have 3 rows and 3 columns. Contact matrices can be used for other groups as long as the dimension of the matrix matches the number of groups.
+
+For example, we might have a meta population model with two geographic areas. Then our contact matrix would be a 2 x 2 matrix with entries representing the contact between and within the geographic areas.
+
+
+## Summary 
+
+In this tutorial, we have learnt the definition of the contact matrix, how they are estimated and how to access social contact data from `socialmixr`. In the next tutorial, we will learn how to use the R package `{epidemics}` to generate disease trajectories from mathematical models with contact matrices from `socialmixr`.
+
+::::::::::::::::::::::::::::::::::::: keypoints 
+
+- `socialmixr` can be used to estimate contact matrices from survey data
+- Contact matrices can be used in different types of analyses
+
+::::::::::::::::::::::::::::::::::::::::::::::::

episodes/modelling-interventions.Rmd

@@ -357,8 +357,8 @@ The equations describing this model are as follows:
 
 $$
 \begin{aligned}
-\frac{dS_i}{dt} & = - \beta S_i \sum_j C_{i,j} I_j -\nu_{t} S_i \\
-\frac{dE_i}{dt} &= \beta S_i\sum_j C_{i,j} I_j - \alpha E_i \\
+\frac{dS_i}{dt} & = - \beta S_i \sum_j C_{i,j} I_j/N_j -\nu_{t} S_i \\
+\frac{dE_i}{dt} &= \beta S_i\sum_j C_{i,j} I_j/N_j - \alpha E_i \\
 \frac{dI_i}{dt} &= \alpha E_i - \gamma I_i \\
 \frac{dR_i}{dt} &=\gamma I_i \\
 \frac{dV_i}{dt} & =\nu_{i,t} S_i\\

episodes/simulating-transmission.Rmd

@@ -211,13 +211,13 @@ For each disease state ($S$, $E$, $I$ and $R$) and age group ($i$), we have a di
 
 $$
 \begin{aligned}
-\frac{dS_i}{dt} & = - \beta S_i \sum_j C_{i,j} I_j \\
-\frac{dE_i}{dt} &= \beta S_i\sum_j C_{i,j} I_j - \alpha E_i \\
+\frac{dS_i}{dt} & = - \beta S_i \sum_j C_{i,j} I_j/N_j \\
+\frac{dE_i}{dt} &= \beta S_i\sum_j C_{i,j} I_j/N_j - \alpha E_i \\
 \frac{dI_i}{dt} &= \alpha E_i - \gamma I_i \\
 \frac{dR_i}{dt} &=\gamma I_i \\
 \end{aligned}
 $$
-Individuals in age group ($i$) move from the susceptible state ($S_i$) to the exposed state ($E_i$) via age-specific contacts with infectious individuals in all groups $\beta S_i \sum_j C_{i,j} I_j$. The contact matrix $C$ allows for heterogeneity in contacts between age groups. They then move to the infectious state at a rate $\alpha$ and recover at a rate $\gamma$. There is no loss of immunity (there are no flows out of the recovered state).
+Individuals in age group ($i$) move from the susceptible state ($S_i$) to the exposed state ($E_i$) via age-specific contacts with infectious individuals in all groups $\beta S_i \sum_j C_{i,j} I_j/N_j$. The contact matrix $C$ allows for heterogeneity in contacts between age groups. They then move to the infectious state at a rate $\alpha$ and recover at a rate $\gamma$. There is no loss of immunity (there are no flows out of the recovered state).
 
 The model parameters are:
 
@@ -248,7 +248,7 @@ To generate trajectories using our model, we must prepare the following inputs:
 
 ### 1. Contact matrix
 
-Contact matrices can be estimated from surveys or contact data, or synthetic ones can be used. We will use the R package `{socialmixr}` to load a contact matrix estimated from POLYMOD survey data [(Mossong et al. 2008)](https://doi.org/10.1371/journal.pmed.0050074).
+We will use the R package `{socialmixr}` to load a contact matrix estimated from POLYMOD survey data [(Mossong et al. 2008)](https://doi.org/10.1371/journal.pmed.0050074).
 
 
 ::::::::::::::::::::::::::::::::::::: challenge 
@@ -297,13 +297,6 @@ contact_matrix
 
 The result is a square matrix with rows and columns for each age group. Contact matrices can be loaded from other sources, but they must be formatted as a matrix to be used in `epidemics`.
 
-::::::::::::::::::::::::::::::::::::: callout
-### Why would a contact matrix be non-symmetric?
-
-One of the arguments of the function `contact_matrix()` is `symmetric=TRUE`. This means that the total number of contacts of age group 1 with age group 2, should be the same as the total number of contacts of age group 2 and age group 1 (see the `socialmixr` [vignette](https://cran.r-project.org/web/packages/socialmixr/vignettes/socialmixr.html) for more detail). However, when contact matrices are estimated from surveys or other sources, the *reported* number of contacts may differ by age group resulting in a non-symmetric contact matrix [(Prem et al 2021)](https://doi.org/10.1371/journal.pcbi.1009098).
-::::::::::::::::::::::::::::::::::::::::::::::::
-
-
 ### 2. Initial conditions
 
 The initial conditions are the proportion of individuals in each disease state $S$, $E$, $I$ and $R$ for each age group at time 0. In this example, we have three age groups age between 0 and 20 years, age between 20 and 40 years and over. Let's assume that in the youngest age category, one in a million individuals are infectious, and the remaining age categories are infection free.