describe the expected qualitative effect of the intervention on epidemic dynamics.
5.2 Setup
-
We will use the diagram package to draw box-and-arrow compartment diagrams, as in Section 6.6.1.2, and the purrr package to create uniform compartment boxes.
+
We will use the diagram package to draw box-and-arrow compartment diagrams, as in R Session 02, and the purrr package to create uniform compartment boxes.
Code
library(diagram)library(purrr)
@@ -583,18 +587,64 @@
\(R\): recovered or removed individuals.
For the standard SIR model, infection moves people from \(S\) to \(I\) at rate \(\lambda\), and recovery or removal moves people from \(I\) to \(R\) at rate \(\gamma\).
-
-5.4 Group activity
-
Choose one intervention or modeling feature that your group wants to represent. Then change the SIR diagram to include it.
-
Examples include:
+
+5.3.1 Exercise 1: Choose one intervention or modeling feature that your group wants to represent. Then change the SIR diagram to include it.
+
Examples include:
vaccination;
isolation or quarantine;
@@ -612,87 +662,137 @@
what parameter or rate each new arrow represents;
what assumptions the new structure makes; and
what qualitative effect you expect the intervention to have on the epidemic curve.
-
-5.5 Example: adding vaccination
+
+
+5.3.2 Example: adding vaccination
One simple way to represent vaccination is to add a flow from \(S\) to \(R\). This assumes vaccination gives protection similar to recovery or removal. That is a strong assumption, but it is a useful starting point.
In this example, \(v\) is the vaccination rate. Before writing equations, ask what this diagram assumes. For example, does vaccination work immediately? Does everyone have the same access to vaccination? Is vaccine protection perfect? Does protection wane?
-
-5.6 Template for your diagram
+
+5.4 Template for your diagram
You can copy and modify this code chunk to make your own diagram. Add compartments to elpos, add arrows to fromto, and label any new rates.
-Code
# Define the locations of each compartment.
-# Each row is one compartment, and the two numbers give its x- and
-# y-position before scaling.
-elpos <-rbind(
-S =c(1, 1),
-I =c(2, 1),
-R =c(3, 1)
-)
-
-# Rescale the compartment positions so they fit inside the plotting area.
-elpos[, 1] <- (2* elpos[, 1] -1) /6
-elpos[, 2] <-0.5
-
-# Define the arrows between compartments.
-# The numbers refer to row positions in elpos, so SI = c(1, 2) draws an
-# arrow from the first row, S, to the second row, I.
-fromto <-rbind(
-SI =c(1, 2),
-IR =c(2, 3)
-)
-
-# Set small plot margins and open a blank plotting area for the diagram.
-op <-par(mar =c(1, 1, 1, 1))
-diagram::openplotmat(asp =0.35)
-
-# Draw one arrow for each row in fromto.
-for (i inseq_len(nrow(fromto))) {
- diagram::straightarrow(
-to = elpos[fromto[i, 2], ],
-from = elpos[fromto[i, 1], ],
-lwd =2,
-arr.pos =0.65,
-arr.length =0.5
- )
-}
-
-# Draw a labeled box for each compartment.
-purrr::walk(
-rownames(elpos),
-.f =function(.x) {
- diagram::textrect(
- elpos[.x, ],
-0.08,
-0.10,
-lab = .x,
-box.col =gray(0.7),
-shadow.col =gray(0.4),
-shadow.size =0.01,
-cex =2
- )
- }
-)
-
-# Label the arrows with the rates that move individuals between states.
-# Adjust the x- and y-positions if you add new compartments or arrows.
-text(mean(elpos[c("S", "I"), 1]), 0.62, expression(lambda), cex =1.8)
-text(mean(elpos[c("I", "R"), 1]), 0.62, expression(gamma), cex =1.8)
-
-# Restore the previous plotting settings.
-par(op)
+Code
# Define the locations of each compartment.
+# Each row is one compartment, and the two numbers give its x- and
+# y-position before scaling.
+elpos <-rbind(
+S =c(1, 1),
+I =c(2, 1),
+R =c(3, 1)
+)
+
+# Rescale the compartment positions so they fit inside the plotting area.
+elpos[, 1] <- (2* elpos[, 1] -1) /6
+elpos[, 2] <-0.5
+
+# Define the arrows between compartments.
+# The numbers refer to row positions in elpos, so SI = c(1, 2) draws an
+# arrow from the first row, S, to the second row, I.
+fromto <-rbind(
+SI =c(1, 2),
+IR =c(2, 3)
+)
+
+# Set small plot margins and open a blank plotting area for the diagram.
+op <-par(mar =c(1, 1, 1, 1))
+diagram::openplotmat(asp =0.35)
+
+# Draw one arrow for each row in fromto.
+for (i inseq_len(nrow(fromto))) {
+ diagram::straightarrow(
+to = elpos[fromto[i, 2], ],
+from = elpos[fromto[i, 1], ],
+lwd =2,
+arr.pos =0.65,
+arr.length =0.5
+ )
+}
+
+# Draw a labeled box for each compartment.
+purrr::walk(
+rownames(elpos),
+.f =function(.x) {
+ diagram::textrect(
+ elpos[.x, ],
+0.08,
+0.10,
+lab = .x,
+box.col =gray(0.7),
+shadow.col =gray(0.4),
+shadow.size =0.01,
+cex =2
+ )
+ }
+)
+
+# Label the arrows with the rates that move individuals between states.
+# Adjust the x- and y-positions if you add new compartments or arrows.
+text(mean(elpos[c("S", "I"), 1]), 0.62, expression(lambda), cex =1.8)
+text(mean(elpos[c("I", "R"), 1]), 0.62, expression(gamma), cex =1.8)
+
+# Restore the previous plotting settings.
+par(op)
-
-5.7 Deliverable
-
Each group should be ready to share:
+
+5.5 Deliverables
+
By the end of the exercise you should be able to share:
one modified box-and-arrow diagram;
a list of the model states;
@@ -1178,276 +1278,279 @@
Source Code
-
---
-title: "Compartmental Diagram Interactive Session"
-execute:
- echo: true
- warning: false
-editor:
- markdown:
- wrap: sentence
-metadata-files:
- - metadata/matthewferrari.yml
----
-
-In this session, you will start with a simple SIR model diagram and modify it to represent an intervention or transmission process that is interesting to your group.
-The goal is not to make the most complicated model possible.
-Instead, the goal is to make each modeling choice explicit: what compartments are needed, what flows connect them, and what assumptions those flows imply.
-
-## Learning goals
-
-By the end of this session, you should be able to:
-
-1. understand how to translate a public-health intervention into a compartment model structure;
-2. identify the states, flows, and parameters in that structure;
-3. explain which assumptions are introduced by adding or removing arrows; and
-4. describe the expected qualitative effect of the intervention on epidemic dynamics.
-
-## Setup
-
-We will use the `diagram` package to draw box-and-arrow compartment diagrams, as in @sec-age-structure-ngm, and the `purrr` package to create uniform compartment boxes.
-
-```{r}
-library(diagram)
-library(purrr)
-```
-
-## Starting point: a simple SIR diagram
-
-The basic SIR model has three compartments:
-
-- $S$: susceptible individuals;
-- $I$: infectious individuals; and
-- $R$: recovered or removed individuals.
-
-```{r}
-#| echo: false
-#| column: body
-#| out-width: 80%
-
-elpos <-rbind(
-S =c(1, 1),
-I =c(2, 1),
-R =c(3, 1)
-)
-
-elpos[, 1] <- (2* elpos[, 1] -1) /6
-elpos[, 2] <-0.5
-
-fromto <-rbind(
-SI =c(1, 2),
-IR =c(2, 3)
-)
-
-op <-par(mar =c(1, 1, 1, 1))
-diagram::openplotmat(asp =0.35)
-
-for (i inseq_len(nrow(fromto))) {
- diagram::straightarrow(
-to = elpos[fromto[i, 2], ],
-from = elpos[fromto[i, 1], ],
-lwd =2,
-arr.pos =0.65,
-arr.length =0.5
- )
-}
-
-purrr::walk(
-c("S", "I", "R"),
-.f =function(.x) {
- diagram::textrect(
- elpos[.x, ],
-0.08,
-0.10,
-lab = .x,
-box.col =gray(0.7),
-shadow.col =gray(0.4),
-shadow.size =0.01,
-cex =2
- )
- }
-)
-
-text(mean(elpos[c("S", "I"), 1]), 0.62, expression(lambda), cex =1.8)
-text(mean(elpos[c("I", "R"), 1]), 0.62, expression(gamma), cex =1.8)
-
-par(op)
-```
-
-For the standard SIR model, infection moves people from $S$ to $I$ at rate $\lambda$, and recovery or removal moves people from $I$ to $R$ at rate $\gamma$.
-
-## Group activity
-
-Choose one intervention or modeling feature that your group wants to represent.
-Then change the SIR diagram to include it.
-
-Examples include:
-
-- vaccination;
-- isolation or quarantine;
-- masking or behavior change;
-- waning immunity;
-- births and deaths;
-- treatment;
-- age groups or other risk groups; or
-- any intervention relevant to your own research area.
-
-For your modified model, prepare a short explanation of:
-
-1. which compartments you added or removed;
-2. which arrows you added or removed;
-3. what parameter or rate each new arrow represents;
-4. what assumptions the new structure makes; and
-5. what qualitative effect you expect the intervention to have on the epidemic curve.
-
-## Example: adding vaccination
-
-One simple way to represent vaccination is to add a flow from $S$ to $R$.
-This assumes vaccination gives protection similar to recovery or removal.
-That is a strong assumption, but it is a useful starting point.
-
-```{r}
-#| echo: false
-#| column: body
-#| out-width: 80%
-
-elpos <-rbind(
-S =c(1, 2),
-I =c(3, 2),
-R =c(2, 1)
-)
-
-elpos[, 1] <- (2* elpos[, 1] -1) /6
-elpos[, 2] <- (2* elpos[, 2] -1) /4
-
-fromto <-rbind(
-SI =c(1, 2),
-IR =c(2, 3),
-SR =c(1, 3)
-)
-
-op <-par(mar =c(1, 1, 1, 1))
-diagram::openplotmat(asp =0.65)
-
-for (i inseq_len(nrow(fromto))) {
- diagram::straightarrow(
-to = elpos[fromto[i, 2], ],
-from = elpos[fromto[i, 1], ],
-lwd =2,
-arr.pos =0.65,
-arr.length =0.5
- )
-}
-
-purrr::walk(
-c("S", "I", "R"),
-.f =function(.x) {
- diagram::textrect(
- elpos[.x, ],
-0.08,
-0.10,
-lab = .x,
-box.col =gray(0.7),
-shadow.col =gray(0.4),
-shadow.size =0.01,
-cex =2
- )
- }
-)
-
-text(mean(elpos[c("S", "I"), 1]), 0.86, expression(lambda), cex =1.8)
-text(mean(elpos[c("I", "R"), 1]) +0.04, 0.47, expression(gamma), cex =1.8)
-text(mean(elpos[c("S", "R"), 1]) -0.04, 0.47, expression(v), cex =1.8)
-
-par(op)
-```
-
-In this example, $v$ is the vaccination rate.
-Before writing equations, ask what this diagram assumes.
-For example, does vaccination work immediately?
-Does everyone have the same access to vaccination?
-Is vaccine protection perfect?
-Does protection wane?
-
-## Template for your diagram
-
-You can copy and modify this code chunk to make your own diagram.
-Add compartments to `elpos`, add arrows to `fromto`, and label any new rates.
-
-```{r}
-#| eval: false
-
-# Define the locations of each compartment.
-# Each row is one compartment, and the two numbers give its x- and
-# y-position before scaling.
-elpos <-rbind(
-S =c(1, 1),
-I =c(2, 1),
-R =c(3, 1)
-)
-
-# Rescale the compartment positions so they fit inside the plotting area.
-elpos[, 1] <- (2* elpos[, 1] -1) /6
-elpos[, 2] <-0.5
-
-# Define the arrows between compartments.
-# The numbers refer to row positions in elpos, so SI = c(1, 2) draws an
-# arrow from the first row, S, to the second row, I.
-fromto <-rbind(
-SI =c(1, 2),
-IR =c(2, 3)
-)
-
-# Set small plot margins and open a blank plotting area for the diagram.
-op <-par(mar =c(1, 1, 1, 1))
-diagram::openplotmat(asp =0.35)
-
-# Draw one arrow for each row in fromto.
-for (i inseq_len(nrow(fromto))) {
- diagram::straightarrow(
-to = elpos[fromto[i, 2], ],
-from = elpos[fromto[i, 1], ],
-lwd =2,
-arr.pos =0.65,
-arr.length =0.5
- )
-}
-
-# Draw a labeled box for each compartment.
-purrr::walk(
-rownames(elpos),
-.f =function(.x) {
- diagram::textrect(
- elpos[.x, ],
-0.08,
-0.10,
-lab = .x,
-box.col =gray(0.7),
-shadow.col =gray(0.4),
-shadow.size =0.01,
-cex =2
- )
- }
-)
-
-# Label the arrows with the rates that move individuals between states.
-# Adjust the x- and y-positions if you add new compartments or arrows.
-text(mean(elpos[c("S", "I"), 1]), 0.62, expression(lambda), cex =1.8)
-text(mean(elpos[c("I", "R"), 1]), 0.62, expression(gamma), cex =1.8)
-
-# Restore the previous plotting settings.
-par(op)
-```
-
-## Deliverable
-
-Each group should be ready to share:
-
-- one modified box-and-arrow diagram;
-- a list of the model states;
-- a list of the model flows and parameters;
-- one or two key assumptions implied by the diagram; and
-- one prediction for how the intervention changes the epidemic trajectory.
+
---
+title: "Compartmental Diagram Interactive Session"
+execute:
+ echo: true
+ warning: false
+editor:
+ markdown:
+ wrap: sentence
+metadata-files:
+ - metadata/matthewferrari.yml
+---
+
+In this session, you will start with a simple SIR model diagram and modify it to represent an intervention or transmission process that is interesting to your group.
+The goal is not to make the most complicated model possible.
+Instead, the goal is to make each modeling choice explicit: what compartments are needed, what flows connect them, and what assumptions those flows imply.
+
+## Learning goals
+
+By the end of this session, you should be able to:
+
+1. understand how to translate a public-health intervention into a compartment model structure;
+2. identify the states, flows, and parameters in that structure;
+3. explain which assumptions are introduced by adding or removing arrows; and
+4. describe the expected qualitative effect of the intervention on epidemic dynamics.
+
+## Setup
+
+We will use the `diagram` package to draw box-and-arrow compartment diagrams, as in [R Session 02](./r-session-02.qmd), and the `purrr` package to create uniform compartment boxes.
+
+```{r}
+library(diagram)
+library(purrr)
+```
+
+## Starting point: a simple SIR diagram
+
+The basic SIR model has three compartments:
+
+- $S$: susceptible individuals;
+- $I$: infectious individuals; and
+- $R$: recovered or removed individuals.
+
+```{r}
+#| echo: true
+#| column: body
+#| out-width: 80%
+#| code-fold: true
+
+elpos <-rbind(
+S =c(1, 1),
+I =c(2, 1),
+R =c(3, 1)
+)
+
+elpos[, 1] <- (2* elpos[, 1] -1) /6
+elpos[, 2] <-0.5
+
+fromto <-rbind(
+SI =c(1, 2),
+IR =c(2, 3)
+)
+
+op <-par(mar =c(1, 1, 1, 1))
+diagram::openplotmat(asp =0.35)
+
+for (i inseq_len(nrow(fromto))) {
+ diagram::straightarrow(
+to = elpos[fromto[i, 2], ],
+from = elpos[fromto[i, 1], ],
+lwd =2,
+arr.pos =0.65,
+arr.length =0.5
+ )
+}
+
+purrr::walk(
+c("S", "I", "R"),
+.f =function(.x) {
+ diagram::textrect(
+ elpos[.x, ],
+0.08,
+0.10,
+lab = .x,
+box.col =gray(0.7),
+shadow.col =gray(0.4),
+shadow.size =0.01,
+cex =2
+ )
+ }
+)
+
+text(mean(elpos[c("S", "I"), 1]), 0.62, expression(lambda), cex =1.8)
+text(mean(elpos[c("I", "R"), 1]), 0.62, expression(gamma), cex =1.8)
+
+par(op)
+```
+
+For the standard SIR model, infection moves people from $S$ to $I$ at rate $\lambda$, and recovery or removal moves people from $I$ to $R$ at rate $\gamma$.
+
+<div class="exercise">
+
+### Exercise 1: Choose one intervention or modeling feature that your group wants to represent. Then change the SIR diagram to include it.
+
+</div>
+
+Examples include:
+
+- vaccination;
+- isolation or quarantine;
+- masking or behavior change;
+- waning immunity;
+- births and deaths;
+- treatment;
+- age groups or other risk groups; or
+- any intervention relevant to your own research area.
+
+For your modified model, prepare a short explanation of:
+
+1. which compartments you added or removed;
+2. which arrows you added or removed;
+3. what parameter or rate each new arrow represents;
+4. what assumptions the new structure makes; and
+5. what qualitative effect you expect the intervention to have on the epidemic curve.
+
+### Example: adding vaccination
+
+One simple way to represent vaccination is to add a flow from $S$ to $R$.
+This assumes vaccination gives protection similar to recovery or removal.
+That is a strong assumption, but it is a useful starting point.
+
+```{r}
+#| echo: true
+#| column: body
+#| out-width: 80%
+#| code-fold: true
+
+elpos <-rbind(
+S =c(1, 2),
+I =c(3, 2),
+R =c(2, 1)
+)
+
+elpos[, 1] <- (2* elpos[, 1] -1) /6
+elpos[, 2] <- (2* elpos[, 2] -1) /4
+
+fromto <-rbind(
+SI =c(1, 2),
+IR =c(2, 3),
+SR =c(1, 3)
+)
+
+op <-par(mar =c(1, 1, 1, 1))
+diagram::openplotmat(asp =0.65)
+
+for (i inseq_len(nrow(fromto))) {
+ diagram::straightarrow(
+to = elpos[fromto[i, 2], ],
+from = elpos[fromto[i, 1], ],
+lwd =2,
+arr.pos =0.65,
+arr.length =0.5
+ )
+}
+
+purrr::walk(
+c("S", "I", "R"),
+.f =function(.x) {
+ diagram::textrect(
+ elpos[.x, ],
+0.08,
+0.10,
+lab = .x,
+box.col =gray(0.7),
+shadow.col =gray(0.4),
+shadow.size =0.01,
+cex =2
+ )
+ }
+)
+
+text(mean(elpos[c("S", "I"), 1]), 0.86, expression(lambda), cex =1.8)
+text(mean(elpos[c("I", "R"), 1]) +0.04, 0.47, expression(gamma), cex =1.8)
+text(mean(elpos[c("S", "R"), 1]) -0.04, 0.47, expression(v), cex =1.8)
+
+par(op)
+```
+
+In this example, $v$ is the vaccination rate.
+Before writing equations, ask what this diagram assumes.
+For example, does vaccination work immediately?
+Does everyone have the same access to vaccination?
+Is vaccine protection perfect?
+Does protection wane?
+
+## Template for your diagram
+
+You can copy and modify this code chunk to make your own diagram.
+Add compartments to `elpos`, add arrows to `fromto`, and label any new rates.
+
+```{r}
+#| eval: false
+
+# Define the locations of each compartment.
+# Each row is one compartment, and the two numbers give its x- and
+# y-position before scaling.
+elpos <-rbind(
+S =c(1, 1),
+I =c(2, 1),
+R =c(3, 1)
+)
+
+# Rescale the compartment positions so they fit inside the plotting area.
+elpos[, 1] <- (2* elpos[, 1] -1) /6
+elpos[, 2] <-0.5
+
+# Define the arrows between compartments.
+# The numbers refer to row positions in elpos, so SI = c(1, 2) draws an
+# arrow from the first row, S, to the second row, I.
+fromto <-rbind(
+SI =c(1, 2),
+IR =c(2, 3)
+)
+
+# Set small plot margins and open a blank plotting area for the diagram.
+op <-par(mar =c(1, 1, 1, 1))
+diagram::openplotmat(asp =0.35)
+
+# Draw one arrow for each row in fromto.
+for (i inseq_len(nrow(fromto))) {
+ diagram::straightarrow(
+to = elpos[fromto[i, 2], ],
+from = elpos[fromto[i, 1], ],
+lwd =2,
+arr.pos =0.65,
+arr.length =0.5
+ )
+}
+
+# Draw a labeled box for each compartment.
+purrr::walk(
+rownames(elpos),
+.f =function(.x) {
+ diagram::textrect(
+ elpos[.x, ],
+0.08,
+0.10,
+lab = .x,
+box.col =gray(0.7),
+shadow.col =gray(0.4),
+shadow.size =0.01,
+cex =2
+ )
+ }
+)
+
+# Label the arrows with the rates that move individuals between states.
+# Adjust the x- and y-positions if you add new compartments or arrows.
+text(mean(elpos[c("S", "I"), 1]), 0.62, expression(lambda), cex =1.8)
+text(mean(elpos[c("I", "R"), 1]), 0.62, expression(gamma), cex =1.8)
+
+# Restore the previous plotting settings.
+par(op)
+```
+
+## Deliverables
+
+By the end of the exercise you should be able to share:
+
+- one modified box-and-arrow diagram;
+- a list of the model states;
+- a list of the model flows and parameters;
+- one or two key assumptions implied by the diagram; and
+- one prediction for how the intervention changes the epidemic trajectory.
This equation shows the important one-to-one relationship between \(R_0\) and the final epidemic size. Plot the relationship between the total epidemic size and \(R_0\) for the complete range of values between 0 and 1.
-
-Show solutions
-
-Code
p_infec <- (seq(0, 1, by =0.001))
-r0_p <- (log(1- p_infec)) / (-p_infec)
-plot(
-x = p_infec,
-y = r0_p,
-main ="Relationship between final proportion infected and R0",
-xlab ="Final proportion infected",
-ylab ="R0"
-)
-
-
-
-
-
-
-
-
+
9.4 Linear Approximation
The next method we introduce takes advantage of the fact that during the early stages of an outbreak, the number of infected individuals is given approximately as \(I(t) \approx I_0 e^{((R_0 - 1)(\gamma + \mu)t)}\). Taking logarithms of both sides, we have \(\ln(I(t)) \approx \ln(I_0) + (R_0 - 1)(\gamma + \mu)t\), showing that the log of the number of infected individuals is approximately linear in time with a slope that reflects both \(R_0\) and the recovery rate.
This suggests that a simple linear regression fit to the first several data points on a log-scale, corrected to account for \(\gamma\) and \(\mu\), provides a rough and ready estimate of \(R_0\). For flu, we can assume \(\mu =0\) because the epidemic occurred over a time period during which natural mortality is negligible. Further, assuming an infectious period of about 2.5 days, we use \(\gamma = (2.5)^{-1} = 0.4\) for the correction. Fitting to the first four data points, we obtain the slope as follows.
@@ -761,148 +757,13 @@
9.4.1 Exercise 2
Our estimate assumes that boys remained infectious during the natural course of infection. The original report on this epidemic indicates that boys found to have symptoms were immediately confined to bed in the infirmary. The report also indicates that only 1 out of 130 adults at the school exhibited any symptoms. It is reasonable, then, to suppose that transmission in each case ceased once he had been admitted to the infirmary. Supposing admission happened within 24 hours of the onset of symptoms. How does this affect our estimate of \(R_0\)? Twelve hours?
-
-Show solutions
-
-Code
#A: If the cases were isolated after 24 hours, then gamma would be 1/1 = 1, and if the cases were isolated after 12 hours, gamma would be 1/0.5 = 2. R0 would be calculated as the beta coefficient over gamma, below:
-
-r0_g1 <-1.094913/1+1
-r0_g1
-
-
[1] 2.094913
-
-Code
r0_g2 <-1.094913/2+1
-r0_g2
-
-
[1] 1.547457
-
-
-
+
9.4.2 Exercise 3
Biweekly data for outbreaks of measles in three communities in Niamey, Niger are provided in the dataframe niamey. Use this method to obtain estimates of \(R_0\) for measles from the first community assuming that the infectious period is approximately two weeks or \(\frac{14}{365} \approx 0.0384\) years.
-
-Show solutions
-
-Code
niamey_ex3 <- niamey
-niamey_ex3[5, 3] <-0#replace a "NA"
-#the command below organizes the data so it can be plotted and analyzed
-niamey_ex3 <-data.frame(
-biweek =rep(seq(1, 16), 3),
-site =c(rep(1, 16), rep(2, 16), rep(3, 16)),
-cases =c(niamey_ex3[, 1], niamey_ex3[, 2], niamey_ex3[, 3])
-) #define "biweeks"
-
-# As the data are reported every two weeks, this corresponds to the 10th observation. Let’s fit a linear model
-
-
-
-Code
# First let's see what the outbreak looks like for the first community on a linear scale
-plot(
- niamey_ex3$biweek[niamey_ex3$site ==1],
- niamey_ex3$cases[niamey_ex3$site ==1],
-type ='p',
-col = niamey_ex3$site,
-xlab ='Biweek',
-ylab ='Cases'
-)
-lines(niamey_ex3$biweek[niamey_ex3$site ==1], niamey_ex3$cases[niamey_ex3$site ==1])
-
-
-
-
-
-
-
-
-Code
# Now let’s try it on a log scale to see until when the outbreak is roughly linear
-plot(
- niamey_ex3$biweek[niamey_ex3$site ==1],
- niamey_ex3$cases[niamey_ex3$site ==1],
-type ='p',
-col = niamey_ex3$site,
-xlab ='Biweek',
-ylab ='Cases',
-log ='y'
-)
-lines(niamey_ex3$biweek[niamey_ex3$site ==1], niamey_ex3$cases[niamey_ex3$site ==1])
-
-
-
-
-
-
-
-
-Code
# here we create a "week" variable to run the analysis on a weekly
-niamey_ex3$week <- niamey_ex3$biweek *2
-# we use the `head` command to take the first N values of a vector. In this case, we're taking the first 10 values of our outcome (cases) and predictor (time) variables.
-model <-lm(
-log(head(niamey_ex3$cases[niamey_ex3$site ==1], 10)) ~
- (head(niamey_ex3$week[niamey_ex3$site ==1], 10))
-)
-summary(model) #summary statistics for fit model
-
-
-Call:
-lm(formula = log(head(niamey_ex3$cases[niamey_ex3$site == 1],
- 10)) ~ (head(niamey_ex3$week[niamey_ex3$site == 1], 10)))
-
-Residuals:
- Min 1Q Median 3Q Max
--0.51796 -0.07364 0.00790 0.10727 0.36805
-
-Coefficients:
- Estimate Std. Error t value
-(Intercept) 2.59679 0.17563 14.79
-head(niamey_ex3$week[niamey_ex3$site == 1], 10) 0.21960 0.01415 15.52
- Pr(>|t|)
-(Intercept) 4.31e-07 ***
-head(niamey_ex3$week[niamey_ex3$site == 1], 10) 2.96e-07 ***
----
-Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
-
-Residual standard error: 0.2571 on 8 degrees of freedom
-Multiple R-squared: 0.9678, Adjusted R-squared: 0.9638
-F-statistic: 240.8 on 1 and 8 DF, p-value: 2.963e-07
-
-
-
-Code
# Now let's get the slope and display it
-slope <-coef(model)[2] #extract slope parameter
-slope #print to screen
As we ran the model by weeks, the \(\gamma\) value is \(2^{-1}\) and \(\hat R_0 = \hat \beta_1 / \gamma +1\) giving \(\hat R_0=0.2196041/0.5+1 \approx 1.44\).
-
+
9.4.3 Exercise 4
A defect with this method is that it uses only a small fraction of the information that might be available, i.e., the first few data points. Indeed, there is nothing in the method that tells one how many data points to use–this is a matter of judgment. Further, there is a tradeoff in that as more and more data points are used the precision of the estimate increases, but this comes at a cost of additional bias. Plot the estimate of \(R_0\) obtained from \(n=3, 4, 5, ...\) data points against the standard error of the slope from the regression analysis to show this tradeoff.
-
-Show solutions
-
Here, we can use a loop and repeat the regression procedure we used above for varying numbers of initial data points in our model.
-
-Code
slope <-NULL
-se <-NULL
-for (i in3:18) {
- model <-lm(
-log(head(niamey_ex3$cases[niamey_ex3$site ==1], i)) ~
- (head(niamey_ex3$week[niamey_ex3$site ==1], i))
- )
- slope <-c(slope, as.numeric(coef(model)[2]))
- se <-c(se, summary(model)$coefficients[4])
-}
-R0 <- slope /0.5+1
-plot(R0, se, ylab ='Standard error')
-
-
-
-
-
-
-
-
+
9.5 Estimating dynamical parameters with least squares
The objective of the previous exercise was to estimate \(R_0\). Knowing \(R_0\) is critical to understanding the dynamics of any epidemic system. It is, however, a composite quantity and is not sufficient to completely describe the epidemic trajectory. For this, we require estimates for all parameters of the model. In this exercise, we introduce a simple approach to model estimation called least squares fitting, sometimes called trajectory matching. The basic idea is that we find the values of the model parameters that minimize the squared differences between model predictions and the observed data. To demonstrate least squares fitting, we consider an outbreak of measles in Niamey, Niger, reported on by (Grais et al. 2006).
@@ -2116,9 +1977,167 @@
9.9.1 Exercise 5
To make things easier, we have assumed the infectious period is known to be 14 days. In terms of years, \(\text{D} = \frac{14}{365} \approx 0.0384\), and the recovery rate is the inverse i.e., \(\gamma = \frac{14}{365}\). Now, modify the code above to estimate \(\gamma\) and \(\beta\) simultaneously.
-
-Show solutions
-
First, add \(\gamma\) as an optimized model parameter.
+
+9.9.2 Exercise 6
+
What happens if one or both of the other unknowns (\(S_0\) and \(I_0\)) is fixed instead of \(\gamma\)?
+
+9.10 Solutions
+
+9.10.1 Exercise 1
+
+Code
p_infec <- (seq(0, 1, by =0.001))
+r0_p <- (log(1- p_infec)) / (-p_infec)
+plot(
+x = p_infec,
+y = r0_p,
+main ="Relationship between final proportion infected and R0",
+xlab ="Final proportion infected",
+ylab ="R0"
+)
+
+
+
+
+
+
+
+
+9.10.2 Exercise 2
+
+Code
#A: If the cases were isolated after 24 hours, then gamma would be 1/1 = 1, and if the cases were isolated after 12 hours, gamma would be 1/0.5 = 2. R0 would be calculated as the beta coefficient over gamma, below:
+
+r0_g1 <-1.094913/1+1
+r0_g1
+
+
[1] 2.094913
+
+Code
r0_g2 <-1.094913/2+1
+r0_g2
+
+
[1] 1.547457
+
+
+
+9.10.3 Exercise 3
+
+Code
niamey[5, 3] <-0#replace a "NA"
+#the command below organizes the data so it can be plotted and analyzed
+niamey <-data.frame(
+biweek =rep(seq(1, 16), 3),
+site =c(rep(1, 16), rep(2, 16), rep(3, 16)),
+cases =c(niamey[, 1], niamey[, 2], niamey[, 3])
+) #define "biweeks"
+
+# As the data are reported every two weeks, this corresponds to the 10th observation. Let’s fit a linear model
+
+
+
+Code
# First let's see what the outbreak looks like for the first community on a linear scale
+plot(
+ niamey$biweek[niamey$site ==1],
+ niamey$cases[niamey$site ==1],
+type ='p',
+col = niamey$site,
+xlab ='Biweek',
+ylab ='Cases'
+)
+lines(niamey$biweek[niamey$site ==1], niamey$cases[niamey$site ==1])
+
+
+
+
+
+
+
+
+Code
# Now let’s try it on a log scale to see until when the outbreak is roughly linear
+plot(
+ niamey$biweek[niamey$site ==1],
+ niamey$cases[niamey$site ==1],
+type ='p',
+col = niamey$site,
+xlab ='Biweek',
+ylab ='Cases',
+log ='y'
+)
+lines(niamey$biweek[niamey$site ==1], niamey$cases[niamey$site ==1])
+
+
+
+
+
+
+
+
+Code
# here we create a "week" variable to run the analysis on a weekly
+niamey$week <- niamey$biweek *2
+# we use the `head` command to take the first N values of a vector. In this case, we're taking the first 10 values of our outcome (cases) and predictor (time) variables.
+model <-lm(
+log(head(niamey$cases[niamey$site ==1], 10)) ~
+ (head(niamey$week[niamey$site ==1], 10))
+)
+summary(model) #summary statistics for fit model
+
+
+Call:
+lm(formula = log(head(niamey$cases[niamey$site == 1], 10)) ~
+ (head(niamey$week[niamey$site == 1], 10)))
+
+Residuals:
+ Min 1Q Median 3Q Max
+-0.51796 -0.07364 0.00790 0.10727 0.36805
+
+Coefficients:
+ Estimate Std. Error t value Pr(>|t|)
+(Intercept) 2.59679 0.17563 14.79 4.31e-07
+head(niamey$week[niamey$site == 1], 10) 0.21960 0.01415 15.52 2.96e-07
+
+(Intercept) ***
+head(niamey$week[niamey$site == 1], 10) ***
+---
+Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
+
+Residual standard error: 0.2571 on 8 degrees of freedom
+Multiple R-squared: 0.9678, Adjusted R-squared: 0.9638
+F-statistic: 240.8 on 1 and 8 DF, p-value: 2.963e-07
+
+
+
+Code
# Now let's get the slope and display it
+slope <-coef(model)[2] #extract slope parameter
+slope #print to screen
As we ran the model by weeks, the \(\gamma\) value is \(2^{-1}\) and \(\hat R_0 = \hat \beta_1 / \gamma +1\) giving \(\hat R_0=0.2196041/0.5+1 \approx 1.44\).
+
+9.10.4 Exercise 4
+
Here, we can use a loop and repeat the regression procedure we used above for varying numbers of initial data points in our model.
+
+Code
slope <-NULL
+se <-NULL
+for (i in3:18) {
+ model <-lm(
+log(head(niamey$cases[niamey$site ==1], i)) ~
+ (head(niamey$week[niamey$site ==1], i))
+ )
+ slope <-c(slope, as.numeric(coef(model)[2]))
+ se <-c(se, summary(model)$coefficients[4])
+}
+R0 <- slope /0.5+1
+plot(R0, se, ylab ='Standard error')
+
+
+
+
+
+
+
+
+9.10.5 Exercise 5
+
First, let’s add in \(\gamma\) estimation into the sse.sir function and create a new sse.sir function called sse.sir.g
What happens if one or both of the other unknowns (\(S_0\) and \(I_0\)) is fixed instead of \(\gamma\)?
-
-Show solutions
-
First we modify the sse_sir_g function so that the fixed initial condition values are supplied separately from the parameters optimized by optim(). Here we fix \(S_0 = 5000\) and/or \(I_0 = 1\), the same values used as starting values in the previous exercise.
\n```\n\n:::\n:::\n\n\n\n::: {.cell .column-body}\n\n```{.r .cell-code}\nR0_mean_age_contacts_df %>%\n select(-sum_cases) %>%\n # Convert to long data frame for facet plotting\n pivot_longer(-R0, names_to = \"metric\", values_to = \"value\") %>%\n ggplot(aes(x = R0, y = value)) +\n geom_line(color = \"slategray4\") +\n geom_point(shape = 21, size = 5, fill = \"slategray4\", alpha = 0.8) +\n facet_wrap(\n ~metric,\n scales = \"free_y\",\n labeller = as_labeller(c(\n mean_age = \"Mean Age of Infection\",\n prev_perc = \"Prevalence (per 100_000) between 15-35 years\"\n ))\n ) +\n labs(\n x = \"R0\",\n y = \"Value\"\n )\n```\n\n::: {.cell-output-display}\n{width=100%}\n:::\n:::\n\n\n### What do real contact networks look like?\n\nThe POLYMOD study [@mossongSocialContactsMixing2008a] was a journal-based look into the contact network in contemporary European society.\nLet's have a look what these data tell us about the contact structure.\n\n\n::: {.cell}\n\n```{.r .cell-code}\nmossong_cont_net <- rio::import(\n \"https://raw.githubusercontent.com/arnold-c/SISMID-Module-02_2023/main/data/mossong-matrix.csv\"\n)\n# mossong_cont_net <- rio::import(here::here(\"data\", \"mossong-matrix.csv\"))\n\nmossong_ages <- unique(mossong_cont_net$contactor)\nmossong_cont_net$contactor <- ordered(\n mossong_cont_net$contactor,\n levels = mossong_ages\n)\n\nmossong_cont_net$contactee <- ordered(\n mossong_cont_net$contactee,\n levels = mossong_ages\n)\n```\n:::\n\n\nSince contacts are symmetric, we'll need to estimate the symmetric contact matrix.\n\n\n::: {.cell}\n\n```{.r .cell-code}\nmossong_mat <- mossong_cont_net %>%\n pivot_wider(\n names_from = contactor,\n values_from = contact.rate\n ) %>%\n select(-contactee) %>%\n as.matrix()\n\nrownames(mossong_mat) <- mossong_ages\n\n# Create a symmetrical contact matrix\nmossong_mat_sym <- (mossong_mat + t(mossong_mat)) / 2\n```\n:::\n\n\nHere we'll use the `filled.contour` function to visualize the contact matrix, to show you an alternative way of visualizing contact matrices.\nNotices that we are using the raw matrix object, not a long dataframe, as previously.\n\n\n::: {.cell .column-body}\n\n```{.r .cell-code}\nfilled.contour(\n ages,\n ages,\n log10(mossong_mat),\n plot.title = title(\n main = \"Log10 of Raw Contact Rate\",\n xlab = \"Age of Contactor\",\n ylab = \"Age of Contactee\"\n )\n)\n```\n\n::: {.cell-output-display}\n{width=100%}\n:::\n:::\n\n\n\n::: {.cell .column-body}\n\n```{.r .cell-code}\nfilled.contour(\n ages,\n ages,\n log10(mossong_mat_sym),\n plot.title = title(\n main = \"Log10 of Symmetrical Contact Rate\",\n xlab = \"Age of Contactor\",\n ylab = \"Age of Contactee\"\n )\n)\n```\n\n::: {.cell-output-display}\n{width=100%}\n:::\n:::\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nmossong_cont_sums <- tibble(\n age = factor(mossong_ages, levels = mossong_ages),\n contactees = rowSums(mossong_mat),\n contactors = colSums(mossong_mat)\n) %>%\n pivot_longer(-age, names_to = \"type\", values_to = \"total_contacts\")\n```\n:::\n\n\n\n::: {.cell .column-body}\n\n```{.r .cell-code}\nggplot(\n mossong_cont_sums,\n aes(\n x = age,\n y = total_contacts,\n color = type,\n fill = type,\n group = type\n )\n) +\n geom_path(linewidth = 1) +\n geom_point(\n position = \"identity\",\n alpha = 0.8,\n shape = 21,\n size = 4\n ) +\n scale_color_manual(\n values = c(\"slategray4\", \"navy\"),\n labels = c(\"Contactees\", \"Contactors\"),\n aesthetics = c(\"color\", \"fill\")\n ) +\n guides(color = \"none\") +\n labs(\n x = \"Age\",\n y = \"Total contacts\",\n fill = \"Type of contact\"\n ) +\n theme(legend.position = \"bottom\")\n```\n\n::: {.cell-output-display}\n{width=100%}\n:::\n:::\n\n\nWhile this matrix tells us how many contacts are made per year by an individual of each age, it doesn't tell us anything about the probability that a contact results in communication of infection.\nLet's assume that each contact has a constant probability $q$ of resulting in a transmission event.\n\n\n::: {.cell}\n\n```{.r .cell-code}\nq <- 3e-5\nmossong_beta_mat <- q * mossong_mat_sym\n```\n:::\n\n\n\n::: {.cell .column-body}\n\n```{.r .cell-code}\nfilled.contour(\n ages,\n ages,\n log10(mossong_beta_mat),\n plot.title = title(\n main = \"WAIFW matrix based on POLYMOD data\",\n xlab = \"Age\",\n ylab = \"Age\"\n )\n)\n```\n\n::: {.cell-output-display}\n{width=100%}\n:::\n:::\n\n\nNow let's simulate the introduction of such a pathogen into a population characterized by this contact structure.\n\n\n::: {.cell}\n\n```{.r .cell-code}\n# Update the parameters with the POLYMOD-based beta matrix\nmossong_params <- multistage_params\nmossong_params[[\"beta_mat\"]] <- mossong_beta_mat\n\n# Solve the model with the updated parameters\nmossong_sol <- deSolve::ode(\n y = demog_yinit_ages,\n times = seq(0, 200, by = 0.5),\n func = multistage_model,\n parms = mossong_params\n)\n\n# Extract the timeseries of infectious individuals\nmossong_infecteds <- mossong_sol[, 1 + iindex]\n\n# Convert infectious individual timeseries to dataframe for plotting\nmossong_infecteds_df <- tibble(\n time = mossong_sol[, 1],\n Juveniles = apply(mossong_infecteds[, juvies], 1, sum),\n Adults = apply(mossong_infecteds[, adults], 1, sum)\n) %>%\n pivot_longer(\n cols = c(Juveniles, Adults),\n names_to = \"age_group\",\n values_to = \"infections\"\n ) %>%\n mutate(\n age_group = factor(age_group, levels = c(\"Juveniles\", \"Adults\"))\n )\n```\n:::\n\n\n\n::: {.cell .column-body}\n\n```{.r .cell-code}\nggplot(\n mossong_infecteds_df,\n aes(x = time, y = infections, color = age_group)\n) +\n geom_line(linewidth = 1.5) +\n scale_color_manual(\n values = age_group_colors\n ) +\n labs(\n x = \"Time\",\n y = \"Number of infections\",\n color = \"Age group\"\n )\n```\n\n::: {.cell-output-display}\n{width=100%}\n:::\n:::\n\n\nAs before, we can also look at the equilibrium seroprevalence\n\n\n::: {.cell}\n\n```{.r .cell-code}\n# Get last time point\nmossong_equil <- drop(tail(mossong_sol, 1))[-1]\n\n# Calculate number of individuals in each age group at end of simulation\nmossong_equil_n <- mossong_equil[sindex] +\n mossong_equil[iindex] +\n mossong_equil[rindex]\n\n# Calculate equilibrium seroprevalence\nmossong_equil_seroprev <- mossong_equil[rindex] / mossong_equil_n\n\n# Convert to dataframe for plotting\nmossong_equil_seroprev_df <- tibble(\n # We can reuse the ages vectors from before as they are the same\n # as the POLYMOD data\n age = ages,\n seroprev = mossong_equil_seroprev,\n width = da_ages\n)\n```\n:::\n\n\n\n::: {.cell .column-body}\n\n```{.r .cell-code}\nggplot(mossong_equil_seroprev_df, aes(x = age, y = seroprev, fill = age)) +\n geom_col(\n width = mossong_equil_seroprev_df$width,\n just = 1.0,\n color = \"black\"\n ) +\n labs(\n x = \"Age\",\n y = \"Seroprevalence\"\n ) +\n scale_x_continuous(breaks = seq(0, 80, 10)) +\n scale_fill_continuous(\n low = age_group_colors[1],\n high = age_group_colors[2]\n )\n```\n\n::: {.cell-output-display}\n{width=100%}\n:::\n:::\n\n\nand compute the $R_0$ for this infection.\n\n\n::: {.cell}\n\n```{.r .cell-code}\nmossong_stable_n <- solve(\n mossong_params[[\"aging_mat\"]],\n -c(mossong_params[[\"births\"]], rep(0, 29))\n)\n\ncalculate_R0(\n beta_mat = mossong_params[[\"beta_mat\"]],\n stable_n_mat = mossong_stable_n,\n aging_mat = mossong_params[[\"aging_mat\"]],\n recovery = mossong_params[[\"recovery\"]]\n)\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n[1] 7.058675\n```\n\n\n:::\n:::\n\n\n::: {.callout-note title=\"QUESTION\"}\nHow does this R0 value compare to the R0 value obtained from @sec-ex-3?\n:::\n\n",
+ "markdown": "---\nsubtitle: \"Heterogeneity and Age Structure in SIR Models\"\nabstract-title: \"\"\nabstract: |\n *Materials adapted from Helen Wearing and Aaron King [@kingAgeStructuredModels2011]*\nexecute:\n warning: false\nmetadata-files:\n - metadata/matthewferrari.yml\n - metadata/mathjax-packages.yml\n---\n\n# R Session 02\n## Load Packages\n\n\n::: {.cell}\n\n```{.r .cell-code}\nlibrary(diagram)\nlibrary(deSolve)\nlibrary(tidyverse)\nlibrary(gt)\nlibrary(rio)\n```\n:::\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\ntheme_set(theme_minimal())\n```\n:::\n\n\n## A Model With 2 Classes\n\nWe'll start with the simplest mechanistic model of two classes we can think of, which has separate classes for two groups $a$ and $b$. These groups could represent different socioeconomic classes, for example.\n\n\n::: {.cell .column-body}\n::: {.cell-output-display}\n{width=100%}\n:::\n:::\n\n\nWhich can be written in equations as,\n$$\n\\begin{aligned}\n \\frac{\\dd{S_a}}{\\dd{t}} &= -\\lambda_a\\,S_a \\phantom{-\\gamma\\,I_b}\\\\\n \\frac{\\dd{S_b}}{\\dd{t}} &= -\\lambda_b\\,S_b \\phantom{-\\gamma\\,I_b}\\\\\n \\frac{\\dd{I_a}}{\\dd{t}} &= \\phantom{-}\\lambda_a\\,S_a -\\gamma\\,I_a\\\\\n \\frac{\\dd{I_b}}{\\dd{t}} &= \\phantom{-}\\lambda_b\\,S_b-\\gamma\\,I_b\\\\\n \\frac{\\dd{R_a}}{\\dd{t}} &= \\phantom{-\\lambda_a\\,S_b}+\\gamma\\,I_a\\\\\n \\frac{\\dd{R_b}}{\\dd{t}} &= \\phantom{-\\lambda_a\\,S_b}+\\gamma\\,I_b\\\\\n \\end{aligned}\n$$\n\nThe $\\lambda$s denote the group-specific force of infections:\n\n$$\n\\begin{aligned}\n \\lambda_a &= \\beta_{aa}\\,I_a+\\beta_{ab}\\,I_b\\\\\n \\lambda_b &= \\beta_{ba}\\,I_a+\\beta_{bb}\\,I_b\n\\end{aligned}\n$$\n\nIn this model, each population can infect each other but the infection moves through the populations separately.\nLet's simulate such a model.\nTo make things concrete, we'll assume that the transmission rates $\\beta$ are greater within groups than between them.\n\n\n::: {.cell}\n\n```{.r .cell-code}\n# Create a named parameter vector that we can index by name in the model\nab_params <- c(\n beta_within = 0.025,\n beta_between = 0.005,\n recovery = 10\n)\n```\n:::\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\n# Here we set up the ODE model that matches the equations above\nab_model <- function(t, x, p, ...) {\n # Unpack the state variables\n Sa <- x[\"Sa\"]\n Sb <- x[\"Sb\"]\n Ia <- x[\"Ia\"]\n Ib <- x[\"Ib\"]\n\n # Unpack the parameters\n beta_within <- p[\"beta_within\"]\n beta_between <- p[\"beta_between\"]\n recovery <- p[\"recovery\"]\n\n # group A force of infection\n lambda_a <- beta_within * Ia + beta_between * Ib\n\n # group B force of infection\n lambda_b <- beta_within * Ib + beta_between * Ia\n\n # The ODEs\n dSadt <- -lambda_a * Sa\n dSbdt <- -lambda_b * Sb\n dIadt <- lambda_a * Sa - recovery * Ia\n dIbdt <- lambda_b * Sb - recovery * Ib\n dRadt <- recovery * Ia\n dRbdt <- recovery * Ib\n\n # Return the derivatives\n list(c(\n dSadt,\n dSbdt,\n dIadt,\n dIbdt,\n dRadt,\n dRbdt\n ))\n}\n```\n:::\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\n# initial conditions\nab_yinit <- c(Sa = 1000, Sb = 2000, Ia = 1, Ib = 1, Ra = 0, Rb = 0)\n\n# Run the ODE solver from the deSolve package\nab_sol <- deSolve::ode(\n y = ab_yinit,\n times = seq(0, 2, by = 0.001),\n func = ab_model,\n parms = ab_params,\n)\n```\n:::\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nab_df <- ab_sol %>%\n # Convert the solution to a tibble for manipulation\n as_tibble() %>%\n # Create and modify columns\n mutate(\n # Convert all columns into type numeric\n across(everything(), as.numeric),\n # Create new columns to track pop sizes in each group\n Na = Sa + Ia + Ra,\n Nb = Sb + Ib + Rb\n ) %>%\n # Go from a wide to long dataframe for ggplot\n pivot_longer(\n cols = -time,\n names_to = c(\"state\", \"group\"),\n names_sep = 1,\n values_to = \"value\"\n ) %>%\n # Clean pivoted columns for ordered plots\n mutate(\n state = factor(state, levels = c(\"S\", \"I\", \"R\", \"N\")),\n group = paste(\"Group\", str_to_upper(group))\n )\n```\n:::\n\n\n\n::: {.cell .column-body}\n\n```{.r .cell-code}\n# Create a vector of colors to be used throughout the ggplots\nSIRcolors <- c(\"#1f77b4\", \"#ff7f0e\", \"#FF3851\", \"#591099\")\n\nggplot(ab_df, aes(x = time, y = value, color = state)) +\n geom_line(linewidth = 1.5) +\n facet_wrap(~group, scales = \"free_y\") +\n scale_color_manual(\n values = SIRcolors,\n labels = c(\"Susceptible\", \"Infected\", \"Recovered\", \"Total\")\n ) +\n labs(\n x = \"Time\",\n y = \"Number of individuals\",\n color = \"State\"\n ) +\n theme(legend.position = \"bottom\")\n```\n\n::: {.cell-output-display}\n{width=100%}\n:::\n:::\n\n\n::: {.callout-question}\nDespite using the same transmission rates, the epidemic in group B is much larger than in group A.\nWhy do you think this is?\n:::\n\nNow let's plot the proportion of individuals in each state for the two groups.\n\n\n::: {.cell}\n\n```{.r .cell-code}\nab_df_props <- ab_df %>%\n # Remove total pop count as we only want the group-specific values\n filter(state != \"N\") %>%\n mutate(\n # Concatenate the state variable and the group letter for each row\n state_group = paste0(state, str_extract_all(group, \"[^Group ]\")),\n # Factor new variable for nicer plotting\n state_group = factor(\n state_group,\n levels = c(\"RA\", \"RB\", \"IA\", \"IB\", \"SA\", \"SB\")\n )\n ) %>%\n # Group by time and state_group so we can calculate the relevant\n # proportions over time\n group_by(time, state_group) %>%\n mutate(\n prop = value / sum(ab_yinit)\n ) %>%\n ungroup()\n```\n:::\n\n\n\n::: {.cell .column-body}\n\n```{.r .cell-code}\n# Create new vectors of colors as using 6: one of each for A and J groups\nScolors <- RColorBrewer::brewer.pal(3, \"Blues\")[c(2, 3)]\nIcolors <- RColorBrewer::brewer.pal(3, \"Oranges\")[c(2, 3)]\nRcolors <- RColorBrewer::brewer.pal(3, \"Greens\")[c(2, 3)]\n\nggplot(ab_df_props, aes(x = time, y = prop, fill = state_group)) +\n geom_area() +\n scale_fill_manual(\n values = c(Scolors, Icolors, Rcolors),\n limits = c(\"SA\", \"SB\", \"IA\", \"IB\", \"RA\", \"RB\"),\n ) +\n labs(\n x = \"Time\",\n y = \"Proportion of individuals\",\n fill = \"State\"\n ) +\n theme(legend.position = \"bottom\")\n```\n\n::: {.cell-output-display}\n{width=100%}\n:::\n:::\n\n\n## A Model With 2 Age Classes\n\nNote that age is a special kind of heterogeneity in an epidemic model because individuals necessarily move from one class (younger) to another class (older) in a directional fashion that is independent of the infection and recovery process.\n\nWe'll start by introducing age into the model above.\nSo now $a$ becomes juveniles and $b$ becomes adults.\nAnd, independent of the disease process, juveniles (of any category) age into adults.\nAdditionally, new juveniles are added through births (always first susceptible) and old individuals are lost to death.\n\n\n::: {.cell .column-body}\n::: {.cell-output-display}\n{width=100%}\n:::\n:::\n\n\nWe can do this very simply using the same ingredients that go into the basic SIR model.\nIn that model, the waiting times in the S and I classes are exponential.\nLet's assume the same thing about the aging process.\nWe'll also add in births into the juvenile susceptible class and deaths from the adult classes.\n\n$$\n \\begin{aligned}\n \\frac{\\dd{S_J}}{\\dd{t}} &= B -\\lambda_J\\,S_J \\phantom{- \\gamma\\,I_A} -\\alpha\\,S_J \\phantom{-\\mu\\,S_A}\\\\\n \\frac{\\dd{S_A}}{\\dd{t}} &= \\phantom{B} - \\lambda_A\\,S_A \\phantom{- \\gamma\\,I_A} +\\alpha\\,S_J -\\mu\\,S_A\\\\\n \\frac{\\dd{I_J}}{\\dd{t}} &= \\phantom{B} +\\lambda_J\\,S_J - \\gamma\\,I_J -\\alpha\\,I_J \\phantom{-\\mu\\,S_A}\\\\\n \\frac{\\dd{I_A}}{\\dd{t}} &= \\phantom{B} +\\lambda_A\\,S_A - \\gamma\\,I_A + \\alpha\\,I_J - \\mu\\,I_A\\\\\n \\frac{\\dd{R_J}}{\\dd{t}} &= \\phantom{B - \\lambda_J\\,S_A} + \\gamma\\,I_J - \\alpha\\,R_J \\phantom{- \\mu\\,S_A}\\\\\n \\frac{\\dd{R_A}}{\\dd{t}} &= \\phantom{B - \\lambda_J\\,S_A} + \\gamma\\,I_A + \\alpha\\,R_J -\\mu\\,R_A\\\\\n \\end{aligned}\n$$\n\nNow, let's simulate this model, under the same assumptions about transmission rates as above.\n\n\n::: {.cell}\n\n```{.r .cell-code}\n# define the parameters for the demographic model\ndemog_params <- c(\n beta_within = 0.004,\n beta_between = 0.002,\n recovery = 10,\n births = 100,\n # Width of age bands in years\n age_band_j = 20,\n age_band_a = 60\n)\n```\n:::\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\ndemog_model <- function(t, x, p, ...) {\n # Unpack states\n Sj <- x[\"Sj\"]\n Sa <- x[\"Sa\"]\n Ij <- x[\"Ij\"]\n Ia <- x[\"Ia\"]\n Rj <- x[\"Rj\"]\n Ra <- x[\"Ra\"]\n\n # Unpack parameters from vector\n beta_within <- p[\"beta_within\"]\n beta_between <- p[\"beta_between\"]\n recovery <- p[\"recovery\"]\n births <- p[\"births\"]\n # Calculate rate of aging from each age group\n aging_j <- 1 / p[\"age_band_j\"]\n aging_a <- 1 / p[\"age_band_a\"]\n\n # juv. force of infection\n lambda_j <- beta_within * Ij + beta_between * Ia\n\n # adult. force of infection\n lambda_a <- beta_within * Ia + beta_between * Ij\n\n # Calculate the ODEs\n dSjdt <- births - (lambda_j * Sj) - (aging_j * Sj)\n dSadt <- -(lambda_a * Sa) + (aging_j * Sj) - (aging_a * Sa)\n dIjdt <- (lambda_j * Sj) - (recovery * Ij) - (aging_j * Ij)\n dIadt <- (lambda_a * Sa) - (recovery * Ia) + (aging_j * Ij) - (aging_a * Ia)\n dRjdt <- (recovery * Ij) - (aging_j * Rj)\n dRadt <- (recovery * Ia) + (aging_j * Rj) - (aging_a * Ra)\n\n # Return the ODEs\n list(c(\n dSjdt,\n dSadt,\n dIjdt,\n dIadt,\n dRjdt,\n dRadt\n ))\n}\n```\n:::\n\n\nNote that in this function, $\\mu=$ `aging_a` $=$ `1 / p[\"age_band_a\"]`, i.e., death, is just like another age class.\n\n\n::: {.cell}\n\n```{.r .cell-code}\n# initial conditions\ndemog_yinit <- c(Sj = 2000, Sa = 3000, Ij = 0, Ia = 1, Rj = 0, Ra = 1000)\n\n# Solve the demographic model\ndemog_sol <- deSolve::ode(\n y = demog_yinit,\n times = seq(0, 200, by = 0.1),\n func = demog_model,\n parms = demog_params\n)\n\ndemog_df <- demog_sol %>%\n as_tibble() %>%\n mutate(\n across(everything(), as.numeric),\n Nj = Sj + Ij + Rj,\n Na = Sa + Ia + Ra,\n # Calculate total population as need for proportional area plots\n N = Nj + Na\n ) %>%\n pivot_longer(\n cols = -c(time, N),\n names_to = c(\"state\", \"group\"),\n names_sep = 1,\n values_to = \"value\"\n ) %>%\n mutate(\n state = factor(state, levels = c(\"S\", \"I\", \"R\", \"N\")),\n group = paste(\"Group\", str_to_upper(group))\n )\n```\n:::\n\n\n
\n\n### Exercise 1: Use this code to plot the number of susceptible, infected, and recovered individuals over time\n\n
\n\n\n::: {.cell .column-body}\n\n```{.r .cell-code}\nggplot(demog_df, aes(x = time, y = value, color = state)) +\n geom_line(linewidth = 1.5) +\n facet_wrap(\n ~group,\n nrow = 2,\n scales = \"free_y\",\n labeller = as_labeller(c(\n `Group A` = \"Adults\",\n `Group J` = \"Juveniles\"\n ))\n ) +\n scale_color_manual(\n values = SIRcolors,\n labels = c(\"Susceptible\", \"Infected\", \"Recovered\", \"Total\")\n ) +\n labs(\n x = \"Time\",\n y = \"Number of individuals\",\n color = \"State\"\n ) +\n theme(legend.position = \"bottom\")\n```\n\n::: {.cell-output-display}\n{width=100%}\n:::\n:::\n\n\nNote that now that births are replenishing susceptibles, infection persists. The results of the above are plotted here:\n\nNow let's plot the proportion of individuals in each state for the two groups.\n\n\n::: {.cell}\n\n```{.r .cell-code}\n# Calculate the proportions in each state and group at each time point\ndemog_df_props <- demog_df %>%\n filter(state != \"N\") %>%\n mutate(\n state_group = paste0(state, str_extract_all(group, \"[^Group ]\")),\n state_group = factor(\n state_group,\n levels = c(\"RJ\", \"RA\", \"IJ\", \"IA\", \"SJ\", \"SA\")\n )\n ) %>%\n group_by(time, state_group) %>%\n mutate(\n # Calculate the proportion of the total population, not the group pop\n prop = value / N\n ) %>%\n ungroup()\n```\n:::\n\n\n\n::: {.cell .column-body}\n\n```{.r .cell-code}\nggplot(demog_df_props, aes(x = time, y = prop, fill = state_group)) +\n geom_area() +\n scale_fill_manual(\n values = c(Scolors, Icolors, Rcolors),\n limits = c(\"SJ\", \"SA\", \"IJ\", \"IA\", \"RJ\", \"RA\")\n ) +\n labs(\n x = \"Time\",\n y = \"Proportion of individuals\",\n fill = \"State\"\n ) +\n theme(legend.position = \"bottom\")\n```\n\n::: {.cell-output-display}\n{width=100%}\n:::\n:::\n\n\nNow let's plot the equilibrium seroprevalence for each age group.\n\n\n::: {.cell}\n\n```{.r .cell-code}\n# Select the last row (time point) of the data frame\ndemog_equil_seroprev <- tail(demog_df) %>%\n mutate(\n # Calculate the proportion of individuals in each state and age group\n prop = value / sum(value),\n # Relabel groups for plots\n group = case_when(group == \"Group J\" ~ \"Juveniles\", TRUE ~ \"Adults\"),\n group = factor(group, levels = c(\"Juveniles\", \"Adults\")),\n .by = group\n ) %>%\n filter(state == \"R\")\n```\n:::\n\n\n\n::: {.cell .column-body}\n\n```{.r .cell-code}\n# Create vector of colors to distinguish between age groups\nage_group_colors <- c(\"#2980B9\", \"#154360\")\n\nggplot(demog_equil_seroprev, aes(x = group, y = prop, fill = group)) +\n geom_col(position = \"identity\") +\n scale_fill_manual(\n values = age_group_colors\n ) +\n labs(\n x = \"Age group\",\n y = \"Equilibrium seroprevalence\",\n fill = \"Age group\"\n ) +\n theme(legend.position = \"none\")\n```\n\n::: {.cell-output-display}\n{width=100%}\n:::\n:::\n\n\nOne thing we are often interested in is the $R_0$ of a system.\nThe details are beyond the scope of this workshop and are not required to complete the exercises in this worksheet, but we have outlined them in @sec-simple-ngm, particularly in @eq-simple-ngm, at the end of this page.\n\n\n::: {.cell}\n\n:::\n\n\nIn our system, $R_0 =$ 2.66.\n\n## Getting more realistic: adding more age classes\n\nIn the models above, the aging process follows an exponential distribution, which means that whether an individual is 1\\~year old or 10 years old, the chance of them becoming an adult is the same!\nTo improve on this, we can assume that the time a juvenile must wait before becoming an adult follows a gamma distribution.\nThis is equivalent to saying that the waiting time is a sum of some number of exponential distributions.\nThis suggests that we can achieve such a distribution by adding age classes to the model, so that becoming an adult means passing through some number of stages.\nWe'll use 30 age classes, and since they don't have to be of equal duration, we'll assume that they're not.\nSpecifically, we'll have 20 1-yr age classes to take us up to adulthood and break adults into 10 age classes of 5\\~yr duration each. The last age class covers age 66-80.\n\nNow, when we had just two age classes, we could write out each of the equations easily enough, but now that we're going to have 30, we'll need to be more systematic.\nIn particular, we'll need to think of $\\beta$ as a matrix of transmission rates.\nLet's see how to define such a matrix in `R`.\nSo that we don't change too many things all at once, let's keep the same contact structure as in the juvenile-adult model.\n\n\n::: {.cell}\n\n```{.r .cell-code}\n# Set up the parameters for model that incorporates a more realistic age matrix\nages_params <- c(\n beta_j = 0.02,\n beta_a = 0.01,\n beta_aj = 0.01 / 2,\n recovery = 10,\n births = 100\n)\n\n# Create a vector of ages\nages <- c(seq(1, 20, by = 1), seq(25, 65, by = 5), 80)\n\n# Calculate the widths of the age bands\nda_ages <- diff(c(0, ages))\n\n# set up a matrix of contact rates between classes: more contact\n# within juveniles and adults than between\nages_beta_mat <- matrix(nrow = 30, ncol = 30)\n\n# transmission rate for juveniles\nages_beta_mat[1:20, 1:20] <- ages_params[\"beta_j\"]\n\n# transmission rate for adults\nages_beta_mat[21:30, 21:30] <- ages_params[\"beta_a\"]\n\n# lower transmission rate between juveniles and adults\nages_beta_mat[1:20, 21:30] <- ages_params[\"beta_aj\"]\n\n# lower transmission rate between juveniles and adults\nages_beta_mat[21:30, 1:20] <- ages_params[\"beta_aj\"]\n```\n:::\n\n\n\n\n\n::: {.cell .column-body}\n\n```{.r .cell-code}\n# Convert matrix to data frame\ntile_df <- expand.grid(x = 1:30, y = 1:30)\ntile_df$value <- as.vector(ages_beta_mat)\n\n# Convert to factor to treat as discrete categories and define colors\ntile_df$value <- factor(tile_df$value)\nbeta_colors <- c(\"0.005\" = \"#fcae91\", \"0.01\" = \"#de2d26\", \"0.02\" = \"#a50f15\")\n\n# Create tile plot with 3 betas\nggplot(tile_df, aes(x = x, y = y, fill = value)) +\n geom_tile(color = \"white\") +\n scale_fill_manual(values = beta_colors, name = expression(beta)) +\n labs(x = \"Age of Contactor\", y = \"Age of Contactee\") +\n scale_x_continuous(breaks = 1:30, labels = as.character(ages)) +\n scale_y_continuous(breaks = 1:30, labels = as.character(ages))\n```\n\n::: {.cell-output-display}\n{width=100%}\n:::\n:::\n\n\n\nWe'll assume that, at the time of introduction, all children are susceptible, as are adults over 45, but that individuals aged 20--45 have seen the pathogen before and are immune.\nThe vector `yinit` expresses these initial conditions.\n\n\n::: {.cell}\n\n```{.r .cell-code}\n# Create a long vector of initial states with only one\n# initial infection in age 50\ndemog_yinit_ages <- c(\n S = c(rep(100, 20), rep(0, 5), rep(200, 5)),\n I = c(rep(0, 25), 1, rep(0, 4)),\n R = c(rep(0, 20), rep(1000, 5), rep(0, 5))\n)\n```\n:::\n\n\nNote that we're starting out with 1 infected individual in the 26th age class (age 50).\n\nThe codes that follow will be a bit easier to follow if we introduce some indexes that will allow us to pick out certain bits of the `yinit` vector.\n\n\n::: {.cell}\n\n```{.r .cell-code}\n# Create vectors of indices relating to each state\n# Note that there are are 30 age classes, 1-20 are 1 year wide, 21-30 are 5 years wide, taking us up to age\n# (ordered S1-80, I1-80, R1-80)\nsindex <- 1:30\niindex <- 31:60\nrindex <- 61:90\n# Create vectors of indices relating to age group\njuvies <- 1:20\nadults <- 21:30\n```\n:::\n\n\nNow, to capture the aging process, it's convenient to define another matrix to hold the rates of movement between age classes.\nGenerally, this matrix would look like this:\n\n$$\n\\begin{pmatrix}\n -\\alpha_1 & 0 & 0 & \\cdots & 0\\\\\n \\alpha_1 & -\\alpha_2 & 0 & \\cdots & 0\\\\\n 0 & \\alpha_2 & -\\alpha_3 & \\cdots & 0\\\\\n \\vdots & & \\ddots & \\ddots & \\vdots \\\\\n 0 & \\cdots & & \\alpha_{29} & -\\alpha_{30}\\\\\n\\end{pmatrix}\n$${#eq-aging-mat}\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\n# Create a diagonal matrix that holds the rates of aging out of each age class\n# The rows represent the age class you're in, the columns represent the age\n# class you're moving to\naging_mat <- diag(-1 / da_ages)\n\n# Fill in the rates of aging into each age class\naging_mat[row(aging_mat) - col(aging_mat) == 1] <- 1 / head(da_ages, -1)\n```\n:::\n\n\nHave a look at the aging matrix, for example by doing:\n\n\n::: {.cell}\n\n```{.r .cell-code}\n# Move fast through the 1-year age classes - negatives are moves out, positives\n# are moves in. Cannot move between non-adjacent age classes\naging_mat[1:5, 1:5]\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n [,1] [,2] [,3] [,4] [,5]\n[1,] -1 0 0 0 0\n[2,] 1 -1 0 0 0\n[3,] 0 1 -1 0 0\n[4,] 0 0 1 -1 0\n[5,] 0 0 0 1 -1\n```\n\n\n:::\n:::\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\n# Move slowly between the wider age classes\naging_mat[25:30, 25:30]\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n [,1] [,2] [,3] [,4] [,5] [,6]\n[1,] -0.2 0.0 0.0 0.0 0.0 0.00000000\n[2,] 0.2 -0.2 0.0 0.0 0.0 0.00000000\n[3,] 0.0 0.2 -0.2 0.0 0.0 0.00000000\n[4,] 0.0 0.0 0.2 -0.2 0.0 0.00000000\n[5,] 0.0 0.0 0.0 0.2 -0.2 0.00000000\n[6,] 0.0 0.0 0.0 0.0 0.2 -0.06666667\n```\n\n\n:::\n:::\n\n\n\n::: {.cell}\n\n:::\n\n\n\n::: {.cell .column-body}\n\n```{.r .cell-code}\naging_mat %>%\n as.data.frame.table() %>%\n mutate(\n age_recipient = rep(ages, 30),\n # Repeat each age in ages vector 30 times before moving to next\n age_source = rep(ages, each = 30)\n ) %>%\n ggplot(aes(\n x = as.factor(age_source),\n y = as.factor(age_recipient),\n z = Freq\n )) +\n geom_tile(colour = \"grey\", size = 0.4, aes(fill = Freq)) +\n scale_fill_gradientn(\n colours = c(\"red\", \"white\", \"blue\"),\n breaks = c(-1, -0.2, 0, 0.2, 1),\n labels = c(\"-1\", \"-0.2\", \"0\", \"0.2\", \"1\")\n ) +\n labs(x = \"Source Age Group\", y = \"Recipient Age Group\")\n```\n\n::: {.cell-output-display}\n{width=100%}\n:::\n:::\n\n\n
\n\n### Exercise 2: What can you say about its structure? How are the different age groups in contact with each other?\n\n
\n\nNow we can put the pieces together to write a simulator for the age-structured SIR dynamics.\n\n\n::: {.cell}\n\n```{.r .cell-code}\n# Using a list instead of a vector to hold the parameters, as ages_beta_mat and\n# aging are both matrices, so we want to keep them as matrices, rather than\n# flattening\nmultistage_params <- list(\n beta_mat = ages_beta_mat,\n recovery = ages_params[\"recovery\"],\n births = ages_params[\"births\"],\n aging_mat = aging_mat\n)\n\nmultistage_model <- function(t, x, p, ...) {\n # Unpack all states from the vector using the relevant indices\n s <- x[sindex]\n i <- x[iindex]\n r <- x[rindex]\n\n # Unpack parameters\n beta_mat <- p[[\"beta_mat\"]]\n recovery <- p[[\"recovery\"]]\n births <- p[[\"births\"]]\n aging_mat <- p[[\"aging_mat\"]]\n\n # Calculate force of infection using matrix multiplication\n lambda <- beta_mat %*% i\n\n # Calculate the ODEs at every time step\n # Note that R add element-wise for vectors i.e. lambda * s results\n # in a vector length 30 (30 age groups), as does aging_mat %*% s,\n # so v1[i] + v2[i] for i in 1:30\n dsdt <- -lambda * s + aging_mat %*% s\n didt <- lambda * s + aging_mat %*% i - recovery * i\n drdt <- aging_mat %*% r + recovery * i\n # Add the birth rate to the first age group\n dsdt[1] <- dsdt[1] + births\n\n # Return the ODEs in a list\n list(c(dsdt, didt, drdt))\n}\n```\n:::\n\n\nWe can plug this into `ode` just as we did the simpler models to simulate an epidemic.\nWe'll then plot the epidemic curve.\n\n\n::: {.cell}\n\n```{.r .cell-code}\n# Solve the model with a realistic age matrix\nmultistage_sol <- deSolve::ode(\n y = demog_yinit_ages,\n times = seq(0, 100, by = 0.1),\n func = multistage_model,\n parms = multistage_params\n)\n\n# Extract all infected age groups at all time points into a new vector\nmultistage_infecteds <- multistage_sol[, 1 + iindex]\n```\n:::\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\n# Create a dataframe of the sum of infectious individuals in Juv/Adult age groups\n# at each time point\nmultistage_df <- tibble(\n # Get all times from model run\n time = multistage_sol[, 1],\n # At each timepoint, apply the sum function to all juvenile infected\n # individuals\n Juveniles = apply(multistage_infecteds[, juvies], 1, sum),\n # At each timepoint, apply the sum function to all adult infected\n # individuals\n Adults = apply(multistage_infecteds[, adults], 1, sum)\n) %>%\n # Pivot to create a long dataframe that works with ggplot\n pivot_longer(\n cols = c(Juveniles, Adults),\n names_to = \"age_group\",\n values_to = \"infections\"\n ) %>%\n # Turn new pivoted variable into a factor to plot nicely\n mutate(\n age_group = factor(age_group, levels = c(\"Juveniles\", \"Adults\"))\n )\n```\n:::\n\n\n\n::: {.cell .column-body}\n\n```{.r .cell-code}\nggplot(multistage_df, aes(x = time, y = infections, color = age_group)) +\n geom_line(linewidth = 1.5) +\n scale_color_manual(\n values = age_group_colors\n ) +\n labs(\n x = \"Time\",\n y = \"Number of infections\",\n color = \"Age group\"\n )\n```\n\n::: {.cell-output-display}\n{width=100%}\n:::\n:::\n\n\nLet's mimic a situation where we have cross-sectional seroprevalence data (e.g. measures of antibodies that tell you someone is in the R class).\nIn using such data, we'd typically assume that the system was at equilibrium.\n\n
\n\n### Exercise 3: What does the equilibrium age-specific seroprevalence look like in this example? {#sec-ex-3}\n\n
\n\nUse the code below to display the age-specific seroprevalence (i.e., the seroprevalence for each age group at equilibrium)\n\n\n::: {.cell}\n\n```{.r .cell-code}\n# Get the last values for all individuals. drop() removes the column name,\n# [-1] removes the time value\nmultistage_equil <- drop(tail(multistage_sol, 1))[-1]\n# Calculate the equilibrium pop sizes of each age group\nmultistage_equil_n <- multistage_equil[sindex] +\n multistage_equil[iindex] +\n multistage_equil[rindex]\n\n# Calculate equilibrium seroprevalence for each age group\nmultistage_equil_seroprev <- multistage_equil[rindex] / multistage_equil_n\n\n# Create a dataframe to store equilibrium seroprev for plotting\nmultistage_equil_seroprev_df <- tibble(\n age = ages,\n seroprev = multistage_equil_seroprev,\n width = da_ages\n)\n```\n:::\n\n\n\n::: {.cell .column-body}\n\n```{.r .cell-code}\nggplot(multistage_equil_seroprev_df, aes(x = age, y = seroprev, fill = age)) +\n # Set column width to width of age bands, and justify to start at\n # lower bound\n geom_col(\n width = multistage_equil_seroprev_df$width,\n just = 1.0,\n color = \"black\"\n ) +\n labs(\n x = \"Age\",\n y = \"Seroprevalence\"\n ) +\n scale_x_continuous(breaks = seq(0, 80, 10)) +\n scale_fill_continuous(\n low = age_group_colors[1],\n high = age_group_colors[2]\n )\n```\n\n::: {.cell-output-display}\n{width=100%}\n:::\n:::\n\n\n::: {.callout-note title=\"QUESTION\"}\nAt what age does the seroprevalence reach 75%?\n:::\n\nLet's also compute $R_0$.\nAnd because we've added a lot of age structure, with transitions between the age groups, we can't just copy and paste the previous Next Generation Matrix code (from @sec-simple-ngm).\nAs before, the details of this computation are out of the workshop's scope, but they are outlined in @sec-age-structure-ngm.\nWe have created a function to calculate R0 for an age-structured SIR and have added some comments, but read @sec-age-structure-ngm for the full details and reasoning.\n\n\n::: {.cell}\n\n```{.r .cell-code}\n# Calculate the stable disease-free age distribution.\n# Could also simulate without any infections.\nmultistage_stable_n <- solve(\n aging_mat,\n c(-1 * multistage_params[[\"births\"]], rep(0, 29))\n)\n```\n:::\n\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\n#' Calculate R0\n#'\n#' Calculate the R0 of an SIR model using the next generation matrix approach\n#' 'described in @heffernanPerspectivesBasicReproductive2005\n#'\n#' @param beta_mat A matrix of beta parameter values\n#' @param stable_n_mat A matrix of the stable age distributions\n#' @param aging_mat A matrix of the aging rates between age compartments\n#' @param recovery_rate The recover rate parameter value (of type double)\n#'\n#' @return The R0 value as type double\n#' @examples\n#' calculate_R0(\n#' beta_mat = multistage_params[[\"beta_mat\"]],\n#' stable_n_mat = multistage_stable_n,\n#' aging_mat = multistage_params[[\"aging_mat\"]],\n#' recovery_rate = multistage_params[[\"recovery\"]]\n#')\ncalculate_R0 <- function(beta_mat, stable_n_mat, aging_mat, recovery_rate) {\n # evaluate new inf jac pde at dfe\n f_mat <- beta_mat * stable_n_mat\n\n # set off-diag of non-inf transition jac pde to neg aging of prev age group\n # (use aging matrix as already calculated in correct places)\n v_mat <- -aging_mat\n # Update the diagonal of non-inf transition jac to add recovery rate\n diag(v_mat) <- diag(v_mat) + recovery_rate\n\n ## Alternative method of calculating using age bands directly\n # v_mat <- diag(recovery_rate + 1 / da_ages)\n # v_mat[row(v_mat) - col(v_mat) == 1] <- - 1 / head(da_ages, -1)\n\n # spectral trace\n R0 <- max(Re(eigen(solve(v_mat, f_mat), only.values = TRUE)$values))\n\n return(R0)\n}\n```\n:::\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\ncalculate_R0(\n beta_mat = multistage_params[[\"beta_mat\"]],\n stable_n_mat = multistage_stable_n,\n aging_mat = multistage_params[[\"aging_mat\"]],\n recovery_rate = multistage_params[[\"recovery\"]]\n)\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n[1] 6.991242\n```\n\n\n:::\n:::\n\n\n
\n\n### Exercise 4: Updating the contact matrix\n\n
\n\n::: {.callout-important}\nYou will need to read and edit the following code carefully so that it runs with your updated parameters.\nWe have highlighted the relevant lines in the code chunks, so hopefully you won't miss them, though make sure you do copy all the code!\n:::\n\n#### Change the juvenile-juvenile contact rate to be 0.025.\n\n\n::: {.cell source-line-numbers='2'}\n\n```{.r .cell-code}\nupdate_age_beta_mat <- ages_beta_mat\nupdate_age_beta_mat[1:20, 1:20] <- ?update_age_params <- multistage_params\nupdate_age_params[[\"beta_mat\"]] <- update_age_beta_mat\n```\n:::\n\n\n\n##### Answer: Change the juvenile-juvenile contact rate to be 0.025.\n\n::: {.cell}\n\n```{.r .cell-code code-fold=\"true\"}\nupdate_age_beta_mat <- ages_beta_mat\nupdate_age_beta_mat[1:20, 1:20] <- 0.025\n\nupdate_age_params <- multistage_params\nupdate_age_params[[\"beta_mat\"]] <- update_age_beta_mat\n```\n:::\n\n\n#### Simulate and plot the age-structured SIR dynamics under your assumptions and record how the age-specific seroprevalence has changed.\n\n\n::: {.cell source-line-numbers='5,15-17,21'}\n\n```{.r .cell-code}\nupdate_age_sol <- deSolve::ode(\n y = demog_yinit_ages,\n times = seq(0, 400, by = 0.1),\n func = multistage_model,\n parms = ?\n)\n\n# Get the time series for each infectious age group\nupdate_age_infecteds <- update_age_sol[, 1 + iindex]\n\n# Get the last values in the time series\nupdate_age_equil <- drop(tail(update_age_sol, 1))[-1]\n\n# Calculate the number of individuals in each age group at the final timepoint\nupdate_age_equil_n <- update_age_equil[ ? ] +\n update_age_equil[ ? ] +\n update_age_equil[ ? ]\n\n# Calculate final seroprevalence\n# Hint: You need PREVIOUSLY infected individuals\nupdate_age_equil_seroprev <- update_age_equil[ ? ] / update_age_equil_n\n\nupdate_age_equil_seroprev_df <- tibble(\n age = ages,\n seroprev = update_age_equil_seroprev,\n width = da_ages\n)\n\nggplot(update_age_equil_seroprev_df, aes(x = age, y = seroprev, fill = age)) +\n # Set column width to width of age bands, and justify to start at\n # lower bound\n geom_col(\n width = update_age_equil_seroprev_df$width,\n just = 1.0, color = \"black\"\n ) +\n labs(\n x = \"Age\",\n y = \"Seroprevalence\"\n ) +\n scale_x_continuous(breaks = seq(0, 80, 10)) +\n scale_fill_continuous(\n low = age_group_colors[1],\n high = age_group_colors[2]\n )\n```\n:::\n\n\n##### Answer: Simulate and plot the age-structured SIR dynamics under your assumptions and record how the age-specific seroprevalence has changed.\n\n::: {.cell}\n\n```{.r .cell-code}\nupdate_age_sol <- deSolve::ode(\n y = demog_yinit_ages,\n times = seq(0, 400, by = 0.1),\n func = multistage_model,\n parms = update_age_params\n)\n\n# Get the time series for each infectious age group\nupdate_age_infecteds <- update_age_sol[, 1 + iindex]\n\nupdate_age_equil <- drop(tail(update_age_sol, 1))[-1]\n\nupdate_age_equil_n <- update_age_equil[sindex] +\n update_age_equil[iindex] +\n update_age_equil[rindex]\n\nupdate_age_equil_seroprev <- update_age_equil[rindex] / update_age_equil_n\n\nupdate_age_equil_seroprev_df <- tibble(\n age = ages,\n seroprev = update_age_equil_seroprev,\n width = da_ages\n)\n```\n:::\n\n\n\n\n::: {.callout-note title=\"QUESTION\"}\nAt what age does the seroprevalence reach 75%?\nHow does this compare to the answer in @sec-ex-3?\n:::\n\n#### Compute $R_0$ for your assumptions.\n\nAs described previously, the calculation for $R_0$ is difficult due to all the age categories and transitions.\nUse the `calculate_R0()` function we [defined earlier](#ngm-function) to calculate $R_0$ for our updated system.\n\n\n::: {.cell source-line-numbers='2-5'}\n\n```{.r .cell-code}\ncalculate_R0(\n beta_mat = ?,\n stable_n_mat = ?,\n aging_mat = ?,\n recovery_rate = ?\n)\n```\n:::\n\n\n\n##### Answer: Compute $R_0$ for your assumptions.\n\n\n::: {.cell}\n\n```{.r .cell-code}\nupdate_age_R0 <- round(\n calculate_R0(\n beta_mat = update_age_params[[\"beta_mat\"]],\n stable_n_mat = multistage_stable_n,\n aging_mat = update_age_params[[\"aging_mat\"]],\n recovery_rate = update_age_params[[\"recovery\"]]\n ),\n digits = 2\n)\n```\n:::\n\n\nIf you've done everything correctly, you should get $R_0 =$ 7.29.\nThis is higher than previously.\nDoes this match your intuition, given our changes to the beta matrix?\n\n## R0 and the Mean Age of Infection\n\nTo develop some intuition about the relationship between $R_0$ and the mean age of infection, let's play with an interactive plot.\nWe will assume that the population is completely susceptible and that the force of infection is constant.\nWe'll also assume that there is heterogenous mixing i.e. no age structure.\n\nAs we've seen in [Matt's lecture on age structure](L05_age-structure.qmd), we can calculate the mean age of infection using the equation below:\n\n$$\nA \\approx \\frac{L}{R_E - 1}\n$${#eq-mean-age}\n\nwhere $L$ is the life expectancy $\\left(L = \\frac{1}{\\mu}\\right)$ and $R_E$ is the effective reproductive number ($R_E = R_0 * (1 - p)$ where $p$ is the fraction of individuals vaccinated).\n\n::: {.callout-note collapse=\"false\"}\nSee @sec-mean-age-r-code for code you can run in `R` to investigate the relationship between $R_0$, vaccination coverage, life expectancy, and the mean age of infection.\n:::\n\n
\n\n### Exercise 5: mean age of infection interactions\n\n
\n\n#### When you increase $R_0$ from 2.0 to 4.0, what happens to the mean age of infection?\n##### Is there a linear change? If not, why not?\n#### With $R_0$ to 4.0, approximately what level of vaccination coverage is required for a mean age of infection of 40 years?\n#### Leaving $R_0$ and vaccination coverage the same, decrease the life expectancy to 50 years. What happens to the mean age of infection?\n##### If it changed, why do you think it did?\n\n \n\n```{ojs}\n//| echo: false\ninit_R0 = 2.0\ninit_vacc = 0.0\ninit_lifeexp = 75\n```\n\n```{ojs}\n//| echo: false\nfunction set(input, value) {\n input.value = value;\n input.dispatchEvent(new Event(\"input\", {bubbles: true}));\n}\n```\n\n```{ojs}\n//| echo: false\n//| panel: sidebar\nviewof reset = Inputs.button([\n [\"Reset all sliders\", () => {\n set(viewof R0, init_R0)\n set(viewof vacc, init_vacc)\n set(viewof lifeexp, init_lifeexp)\n }]\n])\nviewof R0 = Inputs.range(\n [1.0, 10.0],\n {value: 2.0, step: 0.01, label: md`${tex`R_0`}`}\n)\n\nviewof vacc = Inputs.range(\n [0.0, 1.0],\n {value: 0.0, step: 0.01, label: \"Vaccination coverage\"}\n)\n\nviewof lifeexp = Inputs.range(\n [50, 100],\n {value: 75, step: 1, label: \"Life expectancy\"}\n)\n\nmd`${tex`R_E = ${Re_str}`}`\nmd`${tex`\\text{Mean age of infection} = ${Re_mean_age_str}`}`\n```\n\n```{ojs}\n//| echo: false\nRe = R0 * (1 - vacc)\nRe_str = Re.toPrecision(4).toLocaleString()\n```\n\n```{ojs}\n//| echo: false\nfunction calc_mean_age(Re, lifeexp) {\n if(Re >= 1) {\n var mean_age = (lifeexp / (Re - 1))\n } else {\n var mean_age = Infinity\n }\n return mean_age\n}\n```\n\n\n```{ojs}\n//| echo: false\nR0_mean_age = calc_mean_age(R0, lifeexp)\nRe_mean_age = calc_mean_age(Re, lifeexp)\nRe_mean_age_str = Re_mean_age.toPrecision(4).toLocaleString()\n```\n\n\n```{ojs}\n//| echo: false\nimport { aq, op } from '@uwdata/arquero'\n```\n\n```{ojs}\n//| echo: false\nfunction calc_mean_age_arr(vacc, lifeexp, R0_min, R0_max, dR0) {\n var R0_sim = R0_min\n\n var R0 = []\n var Re = []\n var R0_mean_age = []\n var Re_mean_age = []\n\n for (R0_sim = R0_min; R0_sim <= R0_max; R0_sim += dR0) {\n var Re_sim = R0_sim * (1 - vacc)\n var R0_mean_age_sim = calc_mean_age(R0_sim, lifeexp)\n var Re_mean_age_sim = calc_mean_age(Re_sim, lifeexp)\n\n R0.push(R0_sim)\n Re.push(Re_sim)\n R0_mean_age.push(R0_mean_age_sim)\n Re_mean_age.push(Re_mean_age_sim)\n }\n\n return {\n Re: aq.table({\n R0: R0,\n mean_age: Re_mean_age\n }).filter((d) => d.mean_age <= 100),\n R0: aq.table({\n R0: R0,\n mean_age: R0_mean_age\n }).filter((d) => d.mean_age <= 100)\n }\n}\n```\n\n```{ojs}\n//| echo: false\nmean_age_arrs = calc_mean_age_arr(vacc, lifeexp, 1.0, 10.0, 0.01)\n```\n\n```{ojs}\n//| echo: false\nmean_age_dots = [({\n arrow_start: R0_mean_age <= 100 ? R0_mean_age : 100,\n arrow_end: Re_mean_age <= 100 ? Re_mean_age : 100,\n R0: R0.toPrecision(3),\n Re: Re.toPrecision(3),\n R0_mean_age,\n Re_mean_age\n})]\n```\n\n```{ojs}\n//| echo: false\n//| panel: fill\n{\n let R0Color = \"#1f77b4\"\n let ReColor = \"#ff7f0e\"\n\n let plot = Plot.plot({\n color: {\n legend: true,\n domain: [\"Unvaccinated\", \"Vaccinated\"],\n range: [\"#1f77b4\", \"#ff7f0e\"]\n },\n style: {fontSize: \"20px\"},\n marginLeft: 65,\n marginTop: 40,\n marginBottom: 55,\n grid: true,\n width: 800,\n height: 670,\n x: {label: \"R0\", domain: [0, 10]},\n y: {label: \"Mean Age of Infection\", domain: [0, 100]},\n marks: [\n Plot.line(mean_age_arrs.Re, {x: \"R0\", y: \"mean_age\", stroke: ReColor, strokeWidth: 6}),\n Plot.line(mean_age_arrs.R0, {x: \"R0\", y: \"mean_age\", stroke: R0Color, strokeWidth: 6}),\n Re_mean_age <= 100 ?\n [\n vacc > 0.00 ? Plot.dot(mean_age_dots, {x: \"R0\", y: \"Re_mean_age\", r: 12, stroke: ReColor, fill: ReColor, fillOpacity: 0.6}) : null,\n Plot.text(\n mean_age_dots,\n {x: \"R0\", y: \"Re_mean_age\", text: (d) => `Re = ${d.Re}`, dx: 55, dy: -25, fontWeight: \"bold\", fill: ReColor}\n )\n ] :\n null,\n R0_mean_age <= 100 ?\n [\n Plot.dot(mean_age_dots, {x: \"R0\", y: \"R0_mean_age\", r: 12, stroke: R0Color, fill: R0Color, fillOpacity: 0.6}),\n Plot.text(\n mean_age_dots,\n {x: \"R0\", y: \"R0_mean_age\", text: (d) => `R0 = ${d.R0}`, dx: -60, dy: 30, fontWeight: \"bold\", fill: R0Color}\n )\n ] :\n null,\n Plot.arrow(mean_age_dots, {x1: \"R0\", x2: \"R0\", y1: \"arrow_start\", y2: \"arrow_end\", strokeWidth: 4, headLength: 5, inset: 15}),\n ]\n });\n\n return plot;\n}\n```\n\n## Bonus Materials\n### Calculating $R_0$ with the Next Generation Matrix\n#### Simple model structure {#sec-simple-ngm}\n\nTo compute $R_0$, we need to know the stable age distribution (the relative proportion in the juvenile and adult age classes) of the population, which we can find by solving for the disease-free equilibrium: $S_J^*=B/\\alpha$ and $S_A^*=B/\\mu$.\nWith the stable age distribution, we can calculate $R_0$ by constructing the next generation matrix.\nThe code below outlines how the next generation matrix is constructed using the $\\alpha$ (aging from juvenile to adult), $\\mu$ (death), $n$ (total births), $\\gamma$ (recovery), $da$ (width of age groups in years), and $\\beta$ (transmission) parameters.\n\n\nThe next generation matrix is a matrix that specifies how many new age-specific infections are generated by a typical infected individual of each age class (in a fully susceptible population).\nFor example, let's consider an infected adult and ask how many new juvenile infections it generates: this is the product of the number of susceptible juveniles (from the stable age distribution), the per capita transmission rate from adults to juveniles and the average duration of infection, i.e. $S_J^* \\times \\beta_{JA} \\times 1/ (\\gamma+\\mu)$.\nThis forms one element of our next generation matrix.\nThe other elements look very similar, except there are extra terms when we consider an infected juvenile because there is a (very small) chance they may age during the infectious period and therefore cause new infections as an adult:\n\n$$\n\\mathrm{NGM} = \\begin{pmatrix}\n \\frac{S_J^* \\beta_{JJ}}{(\\gamma + \\alpha)} +\n \\frac{\\alpha}{(\\gamma+\\mu)} \\frac{S_J^* \\beta_{JA}}{(\\gamma + \\mu)} &\n \\frac{S_J^* \\beta_{JA}}{(\\gamma + \\mu)} \\\\\n \\frac{S_A^* \\beta_{AJ}}{(\\gamma + \\alpha)} +\n \\frac{\\alpha}{(\\gamma + \\mu)} \\frac{S_A^*\\beta_{AA}}{(\\gamma+\\mu)} &\n \\frac{S_A^* \\beta_{AA}}{(\\gamma + \\mu)}\n \\end{pmatrix}\n$$ {#eq-simple-ngm}\n\n$R_0$ can then be computed as the dominant eigenvalue (i.e., the one with the largest real part) of this matrix. Let's take an example from a model with 2 age classes, from above. First, let's define the components of the next generation matrix:\n\n\n::: {.cell}\n\n```{.r .cell-code}\nngm_params <- c(\n beta_within = 0.011,\n beta_between = 0.005,\n age_band_j = 20,\n age_band_a = 60,\n recovery = 10\n)\n\nalpha_ngm <- 1 / ngm_params[\"age_band_j\"]\nmu_ngm <- 1 / ngm_params[\"age_band_a\"]\nn_ngm <- demog_params[\"births\"] / c(alpha_ngm, mu_ngm)\n\nbeta_ngm <- matrix(\n c(\n ngm_params[\"beta_within\"],\n ngm_params[\"beta_between\"],\n ngm_params[\"beta_between\"],\n ngm_params[\"beta_within\"]\n ),\n nrow = 2,\n ncol = 2\n)\n```\n:::\n\n\nThe Next Generation Matrix can be calculated in `R` as:\n\n\n::: {.cell}\n\n```{.r .cell-code}\nngm <- matrix(\n c(\n n_ngm[1] *\n (beta_ngm[1, 1] / (ngm_params[\"recovery\"] + alpha_ngm)) +\n alpha_ngm /\n (ngm_params[\"recovery\"] + mu_ngm) *\n n_ngm[1] *\n beta_ngm[1, 2] /\n (ngm_params[\"recovery\"] + mu_ngm),\n\n n_ngm[2] *\n beta_ngm[2, 1] /\n (ngm_params[\"recovery\"] + alpha_ngm) +\n alpha_ngm /\n (ngm_params[\"recovery\"] + mu_ngm) *\n n_ngm[2] *\n (beta_ngm[2, 2] / (ngm_params[\"recovery\"] + mu_ngm)),\n\n n_ngm[1] * beta_ngm[1, 2] / (ngm_params[\"recovery\"] + mu_ngm),\n\n n_ngm[2] * beta_ngm[2, 2] / (ngm_params[\"recovery\"] + mu_ngm)\n ),\n nrow = 2,\n ncol = 2\n)\n```\n:::\n\n\nWe can then calculate the eigenvalues and eigenvectors of this matrix:\n\n\n::: {.cell}\n\n```{.r .cell-code}\neigen(ngm)\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\neigen() decomposition\n$values\n[1] 7.191869 1.591188\n\n$vectors\n [,1] [,2]\n[1,] -0.1958841 -0.8560336\n[2,] -0.9806271 0.5169202\n```\n\n\n:::\n:::\n\n\nWe can also choose to just output the eigenvalues:\n\n\n::: {.cell}\n\n```{.r .cell-code}\neigen(ngm, only.values = TRUE)\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n$values\n[1] 7.191869 1.591188\n\n$vectors\nNULL\n```\n\n\n:::\n:::\n\n\nFinally, let's print $R_0$:\n\n\n::: {.cell}\n\n```{.r .cell-code}\nmax(\n Re(\n eigen(ngm, only.values = TRUE)$values\n )\n)\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n[1] 7.191869\n```\n\n\n:::\n:::\n\n\n#### Age-structured NGM {#sec-age-structure-ngm}\n\nMany times it would be impractical to write out the NGM: there are often too many compartments in an age-structured model.\nIn this instance, we want to use a slightly different approach, but the underlying principles are the same: each element of the NGM balances the number of new infections expected to be produced with the rates of individuals coming in and out of that compartment.\n\n##### Stable Age Distribution\n\nThe first thing we need, as before, is the stable age distribution i.e., the disease-free equilibrium.\nThere are two ways we can do this:\n\n1. Simulate the model without any infections for a sufficiently long time (simple, but less accurate)\n2. Do the math.\n\n###### Disease-Free Simulation\n\n\n::: {.cell}\n\n```{.r .cell-code}\n# Set up initial conditions without any infections\nmultistage_sonly_yinit <- c(\n S = c(rep(250, 30)),\n I = c(rep(0, 30)),\n R = c(rep(0, 30))\n)\n\n# Solve disease free sim to get dfe\nmultistage_sonly_sol <- deSolve::ode(\n y = multistage_sonly_yinit,\n times = seq(0, 300, by = 1),\n func = multistage_model,\n parms = multistage_params\n)\n\n# Calculate population size at each time point and save to dataframe\nmultistage_sonly_pop <- tibble(\n time = multistage_sonly_sol[, 1],\n pop = apply(multistage_sonly_sol[, -1], 1, sum)\n)\n```\n:::\n\n\n\n::: {.cell .column-body}\n\n```{.r .cell-code}\nggplot(multistage_sonly_pop, aes(x = time, y = pop)) +\n geom_area(fill = SIRcolors[4], alpha = 0.6) +\n labs(\n x = \"Time\",\n y = \"Population size\"\n )\n```\n\n::: {.cell-output-display}\n{width=100%}\n:::\n:::\n\n\n###### Doing the Math\n\nAlternatively, we can get the stable age distribution by finding the population structure that balances the birth, aging, and death processes.\nWe have already seen the aging matrix in @eq-aging-mat, and at equilibrium, we have the matrix equation\n\n$$\n\\begin{pmatrix}\n -\\alpha_1 & 0 & 0 & \\cdots & 0\\\\\n \\alpha_1 & -\\alpha_2 & 0 & \\cdots & 0\\\\\n 0 & \\alpha_2 & -\\alpha_3 & \\cdots & 0\\\\\n \\vdots & & \\ddots & \\ddots & \\vdots \\\\\n 0 & \\cdots & & \\alpha_{29} & -\\alpha_{30}\\\\\n\\end{pmatrix} .\n\\begin{pmatrix}\n n_1 \\\\ n_2 \\\\ n_3 \\\\ \\vdots \\\\ n_{30}\n\\end{pmatrix} +\n\\begin{pmatrix}\n B \\\\ 0 \\\\ 0 \\\\ \\vdots \\\\ 0\n\\end{pmatrix} =\n\\begin{pmatrix}\n 0 \\\\ 0 \\\\ 0 \\\\ \\vdots \\\\ 0\n\\end{pmatrix}\n$$\n\nTo solve this equation in `R`, we can do\n\n\n::: {.cell}\n\n```{.r .cell-code}\n# solve(a, b) solves the equation a %*% x = b for x, so rearrange equation\n# above so b is on the RHS of the equation\nmultistage_stable_n <- solve(\n aging_mat,\n c(-1 * multistage_params[[\"births\"]], rep(0, 29))\n)\n\nmultistage_stable_n\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n [1] 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100\n[16] 100 100 100 100 100 500 500 500 500 500 500 500 500 500 1500\n```\n\n\n:::\n:::\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\n# Check the final pop value of the S-only sim is equal to the sum of the\n# stable age distribution calculated above\nround(tail(multistage_sonly_pop$pop, 1)) == sum(multistage_stable_n)\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n[1] TRUE\n```\n\n\n:::\n:::\n\n\nThe following lines then compute $R_0$ using the next generation matrix method.\nThis calculation comes from a recipe described in detail previously [@diekmann2000mathematical; @heesterbeekBriefHistoryR02002; @bjornstadAdvancedNextGenerationMatrix2018; @heffernanPerspectivesBasicReproductive2005; @hurfordNextgenerationToolsEvolutionary2009] (we would recommend starting with [@bjornstadAdvancedNextGenerationMatrix2018; and @heffernanPerspectivesBasicReproductive2005]).\n\nThe steps below are copied from [@bjornstadAdvancedNextGenerationMatrix2018]\n\n1. Identify all n infected compartments\n2. Construct a n × 1 matrix, $\\mathbf{F}$, that contains expressions for all completely new infections entering each infected compartment\n3. Construct a n × 1 matrix, $\\mathbf{V^−}$, that contains expressions for all losses out of each infected compartment\n4. Construct a n × 1 matrix, $\\mathbf{V^+}$, that contains expressions for all gains into each infected compartment that does not represent new infections but transfers among infectious classes\n5. Construct a n × 1 matrix, $\\mathbf{V} = \\mathbf{V^−} − \\mathbf{V^+}$\n6. Generate two n × n Jacobian matrices $f$ and $v$ that are the partial derivatives of $\\mathbf{F}$ and $\\mathbf{V}$ with respect to the $n$ infectious state variables\n7. Evaluate the matrices at the disease free equilibrium (dfe), and finally\n8. $R_0$ is the spectral trace (greatest non-negative real eigenvalue) of $\\mathbf{fv}^{−1}|_{\\text{dfe}}$.\n\nWorking through these steps looks like this:\n\n1. Our only infected compartments are the $I_i$ states, for each age group ($i \\in [1, 30]$)\nTo start, let's write out our differential equation:\n\n\\begin{equation}\n \\frac{\\dd{I_i}}{\\dd{t}} = \\lambda_i S_i - \\gamma I_i + \\alpha_{i-1} I_{i-1} - \\alpha_i I_i\n\\end{equation}\n\n\n\n2. We'll now calculate $\\mathbf{F}$ and $\\mathbf{f}$\n\n\n\\begin{align*}\n \\mathbf{F} &= \\begin{pmatrix}\n \\lambda_1 S_1 + \\cancelto{0}{\\alpha_0 I_0} \\\\\n \\vdots \\\\\n \\lambda_{30} S_{30} + \\alpha_{29} I_{29}\n \\end{pmatrix} \\\\ \\\\\n \\mathbf{F} &= \\begin{pmatrix}\n \\left(\\beta_{1, 1} I_1 + \\cdots + \\beta_{1, 30} I_{30} \\right)S_1 \\\\\n \\vdots \\\\\n \\left(\\beta_{30, 1} I_1 + \\cdots + \\beta_{30, 30} I_{30} \\right) S_{30} + \\alpha_{29} I_{29}\n \\end{pmatrix} \\\\\n \\mathbf{f} &= \\begin{pmatrix}\n \\frac{\\partial F_1}{\\partial I_1} & \\cdots & \\frac{\\partial F_1}{\\partial I_{30}} \\\\\n \\vdots & \\ddots & \\vdots \\\\\n \\frac{\\partial F_{30}}{\\partial I_1} & \\cdots & \\frac{\\partial F_{30}}{\\partial I_{30}}\n \\end{pmatrix} & \\frac{\\partial F_1}{\\partial I_1} &= \\frac{\\partial}{\\partial I_1} \\left( \\beta_{1, 1} I_1 + \\cancelto{0}{\\beta_{1, 2} I_2 + \\cdots + \\beta_{1, 30} I_{30}}\\right) S_1 \\\\\n & & \\frac{\\partial F_1}{\\partial I_1} &= \\beta_{1, 1} S_1 \\\\\n \\mathbf{f} &= \\begin{pmatrix}\n \\beta_{1, 1} S_1 & \\cdots & \\beta_{1, 30} S_1 \\\\\n \\vdots & \\ddots & \\vdots \\\\\n \\beta_{30, 1} S_{30} & \\cdots & \\beta_{30, 30} S_{30}\n \\end{pmatrix}\n\\end{align*}\n\n3. Now let's calculate $\\mathbf{V^-}$, $\\mathbf{V^+}$, $\\mathbf{V}$, and $\\mathbf{v}$\n\n\\begin{align*}\n \\mathbf{V^-} &= \\begin{pmatrix}\n \\gamma I_1 + \\alpha_1 I_1 \\\\\n \\gamma I_2 + \\alpha_2 I_2 \\\\\n \\vdots \\\\\n \\gamma I_{30} + \\alpha_{30} I_{30}\n \\end{pmatrix} & \\mathbf{V^+} &= \\begin{pmatrix}\n \\cancelto{0}{\\alpha_0 I_0} \\\\\n \\alpha_1 I_1 \\\\\n \\vdots \\\\\n \\alpha_{29} I_{29}\n \\end{pmatrix} \\\\ \\\\\n \\mathbf{V} &= \\mathbf{V^-} - \\mathbf{V^+} = \\begin{pmatrix}\n \\gamma I_1 + \\alpha_1 I_1 \\\\\n \\gamma I_2 + \\alpha_2 I_2 - \\alpha_1 I_1\\\\\n \\vdots \\\\\n \\gamma I_{30} + \\alpha_{30} I_{30} - \\alpha_{29} I_{29}\n \\end{pmatrix} \\\\ \\\\\n \\mathbf{v} &= \\begin{pmatrix}\n \\pdv{V_1}{I_1} & \\cdots & \\pdv{V_1}{I_{30}} \\\\\n \\vdots & \\ddots & \\vdots \\\\\n \\pdv{V_{30}}{I_1} & \\cdots & \\pdv{V_{30}}{I_{30}}\n \\end{pmatrix} & \\pdv{V_1}{I_1} &= \\pdv{I_1} I_1 \\left( \\gamma + \\alpha_1 \\right)\\\\\n & & \\pdv{V_1}{I_1} &= \\gamma + \\alpha_1 \\\\ \\\\\n & & \\pdv{V_2}{I_1} &= \\pdv{I_1} \\left( \\cancelto{0}{I_2 \\left( \\gamma + \\alpha_2 \\right)} - \\alpha_1 I_1 \\right)\\\\\n & & \\pdv{V_2}{I_1} &= - \\alpha_1 \\\\ \\\\\n & & \\pdv{V_1}{I_2} &= \\pdv{I_2} \\cancelto{0}{I_1 \\left( \\gamma + \\alpha_1 \\right)} \\\\\n & & \\pdv{V_1}{I_2} &= 0 \\\\ \\\\\n \\mathbf{v} &= \\begin{pmatrix}\n \\gamma + \\alpha_1 & 0 & \\cdots & 0\\\\\n - \\alpha_1 & \\gamma + \\alpha_2 & \\cdots & 0 \\\\\n \\vdots & \\ddots & \\ddots & \\vdots \\\\\n 0 & \\cdots & - \\alpha_{29} & \\gamma + \\alpha_{30}\n \\end{pmatrix}\n\\end{align*}\n\n4. To evaluate $\\mathbf{f}$ and $\\mathbf{v}$ at the disease-free equilibrium, we can use the results from our previous calculations.\n$\\mathbf{v}$ doesn't have any state terms in the equation, so it is already evaluated at $\\text{dfe}$.\n$\\mathbf{f}|_{\\text{dfe}}$ involves subsituting $S_i$ for the equilibrium population distribution that balances the births and aging processes.\n\nThis translates to the function we [defined earlier](#ngm-function).\n\n### Mean age of infection `R` code {#sec-mean-age-r-code}\n\nNow let's look at how we can investigate our the relationships between the mean age of infection and $R_0$ and the vaccination coverage using `R`.\nUnlike the interactive plot that simply uses @eq-mean-age to calculate the mean age of infection, we will use a more realistic age-structured model.\n\nLet's return to the earlier models with an age-class mixing matrix.\nBut this time, we'll calculate $R_0$, the mean age of infection, and the number of cases that occur in individuals between 15-35 years as we increase the contact rate.\n\nRecall from the rubella and congenital rubella syndrome (CRS) example that the risk of severe disease outcomes depends on the risk of infection in reproductive age women (here we'll use individuals between 15 and 35 years as a proxy; in reality we would want to account for the differential rate of reproduction at different ages, including those above 35 years).\nRecall also that increasing vaccination reduces $R_E$ -- for simplicity here, so we don't have to add vaccination into the code, we'll simply change $R_0$ because we already know that will give us outcomes that are dynamically equivalent to increasing the proportion of children born who are vaccinated.\nWe'll then calculate how the mean age of infection changes, and specifically how the absolute number of cases among individuals between the ages of 15-35 (as a proxy for reproductive age women) changes.\nTo do so, we'll make a loop and evaluate the code for each of 10 decreaing levels of mixing (which will reduce $R_0$ and we can interpret as analogous to the reduction in $R_E$ that would result from increasing vaccination).\n\n::: {.callout-note title=\"Note about `map()`\" collapse=\"false\"}\nAs you may have noticed previously, we often use the `map_*()` series of functions.\nWe'll use that again here (`map_dfr()`).\nThe full reasons are too complicated to get into here, but broadly speaking, the `map_*()` functions provide us guarantees over the output of our loops.\nIf it runs, we know that something didn't get silently skipped, and that out output vector/list/dataframe is the same length as the inputs.\nThe same can not be said for `for()` loops, and the base `apply` functions are more awkward to work with as they don't have a consistent syntax across the family of functions.\n\nTo learn more, read [this section](just-enough-r.qmd#sec-map-functions) of our `R` primer.\n:::\n\n\n::: {.cell}\n\n```{.r .cell-code}\n# Create vector of scalings to reduce R0\nscale_contact <- seq(1, .2, length = 10)\n\n# Create a new transmission matrix\nbeta_low <- 0.007\nbeta_medium <- 0.02\nbeta_high <- 0.03\n\nbeta_mat <- matrix(beta_low, nrow = 30, ncol = 30)\nbeta_mat[1:20, 1:20] <- beta_medium\nbeta_mat[6:16, 6:16] <- beta_high\n\nscaled_params <- multistage_params\n```\n:::\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\n# Create a dataframe where each row relates to a different R0 value\nR0_mean_age_contacts_df <- map_dfr(\n # Apply the function to each item in the vector of R0 scaling factors\n .x = scale_contact,\n .f = function(.x) {\n # Scale contacts\n scaled_beta_mat <- beta_mat * .x\n\n # Set up parameters\n scaled_params[[\"beta_mat\"]] <- scaled_beta_mat\n\n # Solve the model\n sol <- deSolve::ode(\n y = demog_yinit_ages,\n times = seq(0, 400, by = 0.1),\n func = multistage_model,\n parms = scaled_params\n )\n\n # Get stable age distribution\n stable_n <- solve(\n scaled_params[[\"aging_mat\"]],\n -c(scaled_params[[\"births\"]], rep(0, 29))\n )\n\n R0 <- calculate_R0(\n beta_mat = scaled_params[[\"beta_mat\"]],\n stable_n_mat = stable_n,\n aging_mat = scaled_params[[\"aging_mat\"]],\n recovery = scaled_params[[\"recovery\"]]\n )\n\n sol_dims <- dim(sol)\n\n final_age_sizes <- sol[sol_dims[1], 2:sol_dims[2]]\n\n # Get final number of infected individuals for each S, I, and R class\n susceptibles <- final_age_sizes[sindex]\n infecteds <- final_age_sizes[iindex]\n recovereds <- final_age_sizes[rindex]\n\n # Calculate mean age of infection\n mean_age <- sum(ages * infecteds / sum(infecteds))\n\n # Calculate sum of cases between 15-35 years recall, from the figures\n # above, that that this is the equilibrium prevalence of infection in\n # these age classes, or the average number of individuals that are\n # infected at any given time in these age classes. Note that we're not\n # differentiating between individuals who can and cannot get pregnant\n # here. So we're making an implicit assumption that there no\n # difference in the risk of rubella infection in these groups so that\n # if prevalence goes up in one group, it goes up in the other.\n sum_cases <- sum(infecteds[15:23])\n\n total_15_35 <- sum(susceptibles[15:23]) +\n sum_cases +\n sum(recovereds[15:23])\n\n # Calculate the prevalence as a proportion per 100000 population\n prev_perc <- sum_cases * 100000 / total_15_35\n\n # Return a dataframe with the values\n return(tibble(R0, mean_age, sum_cases, prev_perc))\n }\n)\n```\n:::\n\n\nNow we can make a table of the results and plot mean age and the sum of cases between 15-35 years of age as a function of $R_0$.\n\n\n::: {.cell .column-body}\n\n```{.r .cell-code}\n# Create a table from the dataframe\ngt(R0_mean_age_contacts_df) %>%\n fmt_number(\n columns = everything(),\n decimals = 2\n ) %>%\n # Relabel the column headers\n cols_label(\n R0 = md(\"**R0**\"),\n mean_age = md(\"**Mean age of infection**\"),\n sum_cases = md(\"**Total cases between 15-35 years**\"),\n prev_perc = md(\"**Prevalence (per 100_000) between 15-35 years**\")\n ) %>%\n # Apply style to the table with gray alternating rows\n opt_stylize(style = 1, color = 'gray') %>%\n # Increate space between columns\n opt_horizontal_padding(scale = 3) %>%\n cols_align(\"center\")\n```\n\n::: {.cell-output-display}\n\n```{=html}\n
\n\n
\n \n
\n
R0
\n
Mean age of infection
\n
Total cases between 15-35 years
\n
Prevalence (per 100_000) between 15-35 years
\n
\n \n \n
6.85
\n
6.08
\n
0.54
\n
25.91
\n
6.24
\n
6.61
\n
0.66
\n
31.49
\n
5.63
\n
7.28
\n
0.81
\n
38.34
\n
5.02
\n
8.15
\n
0.98
\n
46.74
\n
4.41
\n
9.30
\n
1.20
\n
56.95
\n
3.81
\n
10.87
\n
1.45
\n
69.10
\n
3.20
\n
13.11
\n
1.74
\n
82.80
\n
2.59
\n
16.40
\n
2.01
\n
95.94
\n
1.98
\n
21.50
\n
2.11
\n
100.71
\n
1.37
\n
29.72
\n
1.48
\n
70.44
\n \n \n
\n
\n```\n\n:::\n:::\n\n\n\n::: {.cell .column-body}\n\n```{.r .cell-code}\nR0_mean_age_contacts_df %>%\n select(-sum_cases) %>%\n # Convert to long data frame for facet plotting\n pivot_longer(-R0, names_to = \"metric\", values_to = \"value\") %>%\n ggplot(aes(x = R0, y = value)) +\n geom_line(color = \"slategray4\") +\n geom_point(shape = 21, size = 5, fill = \"slategray4\", alpha = 0.8) +\n facet_wrap(\n ~metric,\n scales = \"free_y\",\n labeller = as_labeller(c(\n mean_age = \"Mean Age of Infection\",\n prev_perc = \"Prevalence (per 100_000) between 15-35 years\"\n ))\n ) +\n labs(\n x = \"R0\",\n y = \"Value\"\n )\n```\n\n::: {.cell-output-display}\n{width=100%}\n:::\n:::\n\n\n### What do real contact networks look like?\n\nThe POLYMOD study [@mossongSocialContactsMixing2008a] was a journal-based look into the contact network in contemporary European society.\nLet's have a look what these data tell us about the contact structure.\n\n\n::: {.cell}\n\n```{.r .cell-code}\nmossong_cont_net <- rio::import(\n \"https://raw.githubusercontent.com/arnold-c/SISMID-Module-02_2023/main/data/mossong-matrix.csv\"\n)\n# mossong_cont_net <- rio::import(here::here(\"data\", \"mossong-matrix.csv\"))\n\nmossong_ages <- unique(mossong_cont_net$contactor)\nmossong_cont_net$contactor <- ordered(\n mossong_cont_net$contactor,\n levels = mossong_ages\n)\n\nmossong_cont_net$contactee <- ordered(\n mossong_cont_net$contactee,\n levels = mossong_ages\n)\n```\n:::\n\n\nSince contacts are symmetric, we'll need to estimate the symmetric contact matrix.\n\n\n::: {.cell}\n\n```{.r .cell-code}\nmossong_mat <- mossong_cont_net %>%\n pivot_wider(\n names_from = contactor,\n values_from = contact.rate\n ) %>%\n select(-contactee) %>%\n as.matrix()\n\nrownames(mossong_mat) <- mossong_ages\n\n# Create a symmetrical contact matrix\nmossong_mat_sym <- (mossong_mat + t(mossong_mat)) / 2\n```\n:::\n\n\nHere we'll use the `filled.contour` function to visualize the contact matrix, to show you an alternative way of visualizing contact matrices.\nNotices that we are using the raw matrix object, not a long dataframe, as previously.\n\n\n::: {.cell .column-body}\n\n```{.r .cell-code}\nfilled.contour(\n ages,\n ages,\n log10(mossong_mat),\n plot.title = title(\n main = \"Log10 of Raw Contact Rate\",\n xlab = \"Age of Contactor\",\n ylab = \"Age of Contactee\"\n )\n)\n```\n\n::: {.cell-output-display}\n{width=100%}\n:::\n:::\n\n\n\n::: {.cell .column-body}\n\n```{.r .cell-code}\nfilled.contour(\n ages,\n ages,\n log10(mossong_mat_sym),\n plot.title = title(\n main = \"Log10 of Symmetrical Contact Rate\",\n xlab = \"Age of Contactor\",\n ylab = \"Age of Contactee\"\n )\n)\n```\n\n::: {.cell-output-display}\n{width=100%}\n:::\n:::\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nmossong_cont_sums <- tibble(\n age = factor(mossong_ages, levels = mossong_ages),\n contactees = rowSums(mossong_mat),\n contactors = colSums(mossong_mat)\n) %>%\n pivot_longer(-age, names_to = \"type\", values_to = \"total_contacts\")\n```\n:::\n\n\n\n::: {.cell .column-body}\n\n```{.r .cell-code}\nggplot(\n mossong_cont_sums,\n aes(\n x = age,\n y = total_contacts,\n color = type,\n fill = type,\n group = type\n )\n) +\n geom_path(linewidth = 1) +\n geom_point(\n position = \"identity\",\n alpha = 0.8,\n shape = 21,\n size = 4\n ) +\n scale_color_manual(\n values = c(\"slategray4\", \"navy\"),\n labels = c(\"Contactees\", \"Contactors\"),\n aesthetics = c(\"color\", \"fill\")\n ) +\n guides(color = \"none\") +\n labs(\n x = \"Age\",\n y = \"Total contacts\",\n fill = \"Type of contact\"\n ) +\n theme(legend.position = \"bottom\")\n```\n\n::: {.cell-output-display}\n{width=100%}\n:::\n:::\n\n\nWhile this matrix tells us how many contacts are made per year by an individual of each age, it doesn't tell us anything about the probability that a contact results in communication of infection.\nLet's assume that each contact has a constant probability $q$ of resulting in a transmission event.\n\n\n::: {.cell}\n\n```{.r .cell-code}\nq <- 3e-5\nmossong_beta_mat <- q * mossong_mat_sym\n```\n:::\n\n\n\n::: {.cell .column-body}\n\n```{.r .cell-code}\nfilled.contour(\n ages,\n ages,\n log10(mossong_beta_mat),\n plot.title = title(\n main = \"WAIFW matrix based on POLYMOD data\",\n xlab = \"Age\",\n ylab = \"Age\"\n )\n)\n```\n\n::: {.cell-output-display}\n{width=100%}\n:::\n:::\n\n\nNow let's simulate the introduction of such a pathogen into a population characterized by this contact structure.\n\n\n::: {.cell}\n\n```{.r .cell-code}\n# Update the parameters with the POLYMOD-based beta matrix\nmossong_params <- multistage_params\nmossong_params[[\"beta_mat\"]] <- mossong_beta_mat\n\n# Solve the model with the updated parameters\nmossong_sol <- deSolve::ode(\n y = demog_yinit_ages,\n times = seq(0, 200, by = 0.5),\n func = multistage_model,\n parms = mossong_params\n)\n\n# Extract the timeseries of infectious individuals\nmossong_infecteds <- mossong_sol[, 1 + iindex]\n\n# Convert infectious individual timeseries to dataframe for plotting\nmossong_infecteds_df <- tibble(\n time = mossong_sol[, 1],\n Juveniles = apply(mossong_infecteds[, juvies], 1, sum),\n Adults = apply(mossong_infecteds[, adults], 1, sum)\n) %>%\n pivot_longer(\n cols = c(Juveniles, Adults),\n names_to = \"age_group\",\n values_to = \"infections\"\n ) %>%\n mutate(\n age_group = factor(age_group, levels = c(\"Juveniles\", \"Adults\"))\n )\n```\n:::\n\n\n\n::: {.cell .column-body}\n\n```{.r .cell-code}\nggplot(\n mossong_infecteds_df,\n aes(x = time, y = infections, color = age_group)\n) +\n geom_line(linewidth = 1.5) +\n scale_color_manual(\n values = age_group_colors\n ) +\n labs(\n x = \"Time\",\n y = \"Number of infections\",\n color = \"Age group\"\n )\n```\n\n::: {.cell-output-display}\n{width=100%}\n:::\n:::\n\n\nAs before, we can also look at the equilibrium seroprevalence\n\n\n::: {.cell}\n\n```{.r .cell-code}\n# Get last time point\nmossong_equil <- drop(tail(mossong_sol, 1))[-1]\n\n# Calculate number of individuals in each age group at end of simulation\nmossong_equil_n <- mossong_equil[sindex] +\n mossong_equil[iindex] +\n mossong_equil[rindex]\n\n# Calculate equilibrium seroprevalence\nmossong_equil_seroprev <- mossong_equil[rindex] / mossong_equil_n\n\n# Convert to dataframe for plotting\nmossong_equil_seroprev_df <- tibble(\n # We can reuse the ages vectors from before as they are the same\n # as the POLYMOD data\n age = ages,\n seroprev = mossong_equil_seroprev,\n width = da_ages\n)\n```\n:::\n\n\n\n::: {.cell .column-body}\n\n```{.r .cell-code}\nggplot(mossong_equil_seroprev_df, aes(x = age, y = seroprev, fill = age)) +\n geom_col(\n width = mossong_equil_seroprev_df$width,\n just = 1.0,\n color = \"black\"\n ) +\n labs(\n x = \"Age\",\n y = \"Seroprevalence\"\n ) +\n scale_x_continuous(breaks = seq(0, 80, 10)) +\n scale_fill_continuous(\n low = age_group_colors[1],\n high = age_group_colors[2]\n )\n```\n\n::: {.cell-output-display}\n{width=100%}\n:::\n:::\n\n\nand compute the $R_0$ for this infection.\n\n\n::: {.cell}\n\n```{.r .cell-code}\nmossong_stable_n <- solve(\n mossong_params[[\"aging_mat\"]],\n -c(mossong_params[[\"births\"]], rep(0, 29))\n)\n\ncalculate_R0(\n beta_mat = mossong_params[[\"beta_mat\"]],\n stable_n_mat = mossong_stable_n,\n aging_mat = mossong_params[[\"aging_mat\"]],\n recovery = mossong_params[[\"recovery\"]]\n)\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n[1] 7.058675\n```\n\n\n:::\n:::\n\n\n::: {.callout-note title=\"QUESTION\"}\nHow does this R0 value compare to the R0 value obtained from @sec-ex-3?\n:::\n\n",
"supporting": [],
"filters": [
"rmarkdown/pagebreak.lua"
diff --git a/_freeze/r-session-02/figure-html/unnamed-chunk-10-1.png b/_freeze/r-session-02/figure-html/unnamed-chunk-10-1.png
index 261f882..67e5a05 100644
Binary files a/_freeze/r-session-02/figure-html/unnamed-chunk-10-1.png and b/_freeze/r-session-02/figure-html/unnamed-chunk-10-1.png differ
diff --git a/_freeze/r-session-02/figure-html/unnamed-chunk-11-1.png b/_freeze/r-session-02/figure-html/unnamed-chunk-11-1.png
index 2ecb6ee..8deb5be 100644
Binary files a/_freeze/r-session-02/figure-html/unnamed-chunk-11-1.png and b/_freeze/r-session-02/figure-html/unnamed-chunk-11-1.png differ
diff --git a/_freeze/r-session-02/figure-html/unnamed-chunk-15-1.png b/_freeze/r-session-02/figure-html/unnamed-chunk-15-1.png
index 7b07080..9d93011 100644
Binary files a/_freeze/r-session-02/figure-html/unnamed-chunk-15-1.png and b/_freeze/r-session-02/figure-html/unnamed-chunk-15-1.png differ
diff --git a/_freeze/r-session-02/figure-html/unnamed-chunk-17-1.png b/_freeze/r-session-02/figure-html/unnamed-chunk-17-1.png
index 2bd3c42..f07648a 100644
Binary files a/_freeze/r-session-02/figure-html/unnamed-chunk-17-1.png and b/_freeze/r-session-02/figure-html/unnamed-chunk-17-1.png differ
diff --git a/_freeze/r-session-02/figure-html/unnamed-chunk-19-1.png b/_freeze/r-session-02/figure-html/unnamed-chunk-19-1.png
index c302009..0a6dbe0 100644
Binary files a/_freeze/r-session-02/figure-html/unnamed-chunk-19-1.png and b/_freeze/r-session-02/figure-html/unnamed-chunk-19-1.png differ
diff --git a/_freeze/r-session-02/figure-html/unnamed-chunk-22-1.png b/_freeze/r-session-02/figure-html/unnamed-chunk-22-1.png
index 8c44fa5..e9b61bb 100644
Binary files a/_freeze/r-session-02/figure-html/unnamed-chunk-22-1.png and b/_freeze/r-session-02/figure-html/unnamed-chunk-22-1.png differ
diff --git a/_freeze/r-session-02/figure-html/unnamed-chunk-29-1.png b/_freeze/r-session-02/figure-html/unnamed-chunk-29-1.png
index 79008fd..507a5c6 100644
Binary files a/_freeze/r-session-02/figure-html/unnamed-chunk-29-1.png and b/_freeze/r-session-02/figure-html/unnamed-chunk-29-1.png differ
diff --git a/_freeze/r-session-02/figure-html/unnamed-chunk-3-1.png b/_freeze/r-session-02/figure-html/unnamed-chunk-3-1.png
index 723c5f3..d867486 100644
Binary files a/_freeze/r-session-02/figure-html/unnamed-chunk-3-1.png and b/_freeze/r-session-02/figure-html/unnamed-chunk-3-1.png differ
diff --git a/_freeze/r-session-02/figure-html/unnamed-chunk-33-1.png b/_freeze/r-session-02/figure-html/unnamed-chunk-33-1.png
index 1ebb2ef..b590e34 100644
Binary files a/_freeze/r-session-02/figure-html/unnamed-chunk-33-1.png and b/_freeze/r-session-02/figure-html/unnamed-chunk-33-1.png differ
diff --git a/_freeze/r-session-02/figure-html/unnamed-chunk-35-1.png b/_freeze/r-session-02/figure-html/unnamed-chunk-35-1.png
index 5501269..294e964 100644
Binary files a/_freeze/r-session-02/figure-html/unnamed-chunk-35-1.png and b/_freeze/r-session-02/figure-html/unnamed-chunk-35-1.png differ
diff --git a/_freeze/r-session-02/figure-html/unnamed-chunk-62-1.png b/_freeze/r-session-02/figure-html/unnamed-chunk-62-1.png
index c3b6cd1..8e617c7 100644
Binary files a/_freeze/r-session-02/figure-html/unnamed-chunk-62-1.png and b/_freeze/r-session-02/figure-html/unnamed-chunk-62-1.png differ
diff --git a/_freeze/r-session-02/figure-html/unnamed-chunk-68-1.png b/_freeze/r-session-02/figure-html/unnamed-chunk-68-1.png
index 8bfa465..af1da25 100644
Binary files a/_freeze/r-session-02/figure-html/unnamed-chunk-68-1.png and b/_freeze/r-session-02/figure-html/unnamed-chunk-68-1.png differ
diff --git a/_freeze/r-session-02/figure-html/unnamed-chunk-71-1.png b/_freeze/r-session-02/figure-html/unnamed-chunk-71-1.png
index 6cfc3e0..82267c0 100644
Binary files a/_freeze/r-session-02/figure-html/unnamed-chunk-71-1.png and b/_freeze/r-session-02/figure-html/unnamed-chunk-71-1.png differ
diff --git a/_freeze/r-session-02/figure-html/unnamed-chunk-72-1.png b/_freeze/r-session-02/figure-html/unnamed-chunk-72-1.png
index 8bdc722..6c093d9 100644
Binary files a/_freeze/r-session-02/figure-html/unnamed-chunk-72-1.png and b/_freeze/r-session-02/figure-html/unnamed-chunk-72-1.png differ
diff --git a/_freeze/r-session-02/figure-html/unnamed-chunk-74-1.png b/_freeze/r-session-02/figure-html/unnamed-chunk-74-1.png
index d37121d..59262e7 100644
Binary files a/_freeze/r-session-02/figure-html/unnamed-chunk-74-1.png and b/_freeze/r-session-02/figure-html/unnamed-chunk-74-1.png differ
diff --git a/_freeze/r-session-02/figure-html/unnamed-chunk-76-1.png b/_freeze/r-session-02/figure-html/unnamed-chunk-76-1.png
index 71f8f41..efb6be6 100644
Binary files a/_freeze/r-session-02/figure-html/unnamed-chunk-76-1.png and b/_freeze/r-session-02/figure-html/unnamed-chunk-76-1.png differ
diff --git a/_freeze/r-session-02/figure-html/unnamed-chunk-78-1.png b/_freeze/r-session-02/figure-html/unnamed-chunk-78-1.png
index 273f4ba..740a6e8 100644
Binary files a/_freeze/r-session-02/figure-html/unnamed-chunk-78-1.png and b/_freeze/r-session-02/figure-html/unnamed-chunk-78-1.png differ
diff --git a/_freeze/r-session-02/figure-html/unnamed-chunk-8-1.png b/_freeze/r-session-02/figure-html/unnamed-chunk-8-1.png
index 3383441..7c926fb 100644
Binary files a/_freeze/r-session-02/figure-html/unnamed-chunk-8-1.png and b/_freeze/r-session-02/figure-html/unnamed-chunk-8-1.png differ
diff --git a/_freeze/r-session-02/figure-html/unnamed-chunk-80-1.png b/_freeze/r-session-02/figure-html/unnamed-chunk-80-1.png
index 0d0dd58..4f5837a 100644
Binary files a/_freeze/r-session-02/figure-html/unnamed-chunk-80-1.png and b/_freeze/r-session-02/figure-html/unnamed-chunk-80-1.png differ
diff --git a/day-1-interactive-session.qmd b/day-1-interactive-session.qmd
index b5d655b..9d04d2d 100644
--- a/day-1-interactive-session.qmd
+++ b/day-1-interactive-session.qmd
@@ -25,7 +25,7 @@ By the end of this session, you should be able to:
## Setup
-We will use the `diagram` package to draw box-and-arrow compartment diagrams, as in @sec-age-structure-ngm, and the `purrr` package to create uniform compartment boxes.
+We will use the `diagram` package to draw box-and-arrow compartment diagrams, as in [R Session 02](./r-session-02.qmd), and the `purrr` package to create uniform compartment boxes.
```{r}
library(diagram)
@@ -41,9 +41,10 @@ The basic SIR model has three compartments:
- $R$: recovered or removed individuals.
```{r}
-#| echo: false
+#| echo: true
#| column: body
#| out-width: 80%
+#| code-fold: true
elpos <- rbind(
S = c(1, 1),
@@ -96,10 +97,11 @@ par(op)
For the standard SIR model, infection moves people from $S$ to $I$ at rate $\lambda$, and recovery or removal moves people from $I$ to $R$ at rate $\gamma$.
-## Group activity
+
-Choose one intervention or modeling feature that your group wants to represent.
-Then change the SIR diagram to include it.
+### Exercise 1: Choose one intervention or modeling feature that your group wants to represent. Then change the SIR diagram to include it.
+
+
Examples include:
@@ -120,16 +122,17 @@ For your modified model, prepare a short explanation of:
4. what assumptions the new structure makes; and
5. what qualitative effect you expect the intervention to have on the epidemic curve.
-## Example: adding vaccination
+### Example: adding vaccination
One simple way to represent vaccination is to add a flow from $S$ to $R$.
This assumes vaccination gives protection similar to recovery or removal.
That is a strong assumption, but it is a useful starting point.
```{r}
-#| echo: false
+#| echo: true
#| column: body
#| out-width: 80%
+#| code-fold: true
elpos <- rbind(
S = c(1, 2),
@@ -259,9 +262,9 @@ text(mean(elpos[c("I", "R"), 1]), 0.62, expression(gamma), cex = 1.8)
par(op)
```
-## Deliverable
+## Deliverables
-Each group should be ready to share:
+By the end of the exercise you should be able to share:
- one modified box-and-arrow diagram;
- a list of the model states;
diff --git a/r-session-02.qmd b/r-session-02.qmd
index 03ad683..41b1d85 100644
--- a/r-session-02.qmd
+++ b/r-session-02.qmd
@@ -836,7 +836,7 @@ aging_mat %>%
z = Freq
)) +
geom_tile(colour = "grey", size = 0.4, aes(fill = Freq)) +
-scale_fill_gradientn(
+ scale_fill_gradientn(
colours = c("red", "white", "blue"),
breaks = c(-1, -0.2, 0, 0.2, 1),
labels = c("-1", "-0.2", "0", "0.2", "1")
diff --git a/r-session-02_files/figure-html/unnamed-chunk-10-1.png b/r-session-02_files/figure-html/unnamed-chunk-10-1.png
index 261f882..67e5a05 100644
Binary files a/r-session-02_files/figure-html/unnamed-chunk-10-1.png and b/r-session-02_files/figure-html/unnamed-chunk-10-1.png differ
diff --git a/r-session-02_files/figure-html/unnamed-chunk-11-1.png b/r-session-02_files/figure-html/unnamed-chunk-11-1.png
index 2ecb6ee..8deb5be 100644
Binary files a/r-session-02_files/figure-html/unnamed-chunk-11-1.png and b/r-session-02_files/figure-html/unnamed-chunk-11-1.png differ
diff --git a/r-session-02_files/figure-html/unnamed-chunk-15-1.png b/r-session-02_files/figure-html/unnamed-chunk-15-1.png
index 7b07080..9d93011 100644
Binary files a/r-session-02_files/figure-html/unnamed-chunk-15-1.png and b/r-session-02_files/figure-html/unnamed-chunk-15-1.png differ
diff --git a/r-session-02_files/figure-html/unnamed-chunk-17-1.png b/r-session-02_files/figure-html/unnamed-chunk-17-1.png
index 2bd3c42..f07648a 100644
Binary files a/r-session-02_files/figure-html/unnamed-chunk-17-1.png and b/r-session-02_files/figure-html/unnamed-chunk-17-1.png differ
diff --git a/r-session-02_files/figure-html/unnamed-chunk-19-1.png b/r-session-02_files/figure-html/unnamed-chunk-19-1.png
index c302009..0a6dbe0 100644
Binary files a/r-session-02_files/figure-html/unnamed-chunk-19-1.png and b/r-session-02_files/figure-html/unnamed-chunk-19-1.png differ
diff --git a/r-session-02_files/figure-html/unnamed-chunk-22-1.png b/r-session-02_files/figure-html/unnamed-chunk-22-1.png
index 8c44fa5..e9b61bb 100644
Binary files a/r-session-02_files/figure-html/unnamed-chunk-22-1.png and b/r-session-02_files/figure-html/unnamed-chunk-22-1.png differ
diff --git a/r-session-02_files/figure-html/unnamed-chunk-29-1.png b/r-session-02_files/figure-html/unnamed-chunk-29-1.png
index 79008fd..507a5c6 100644
Binary files a/r-session-02_files/figure-html/unnamed-chunk-29-1.png and b/r-session-02_files/figure-html/unnamed-chunk-29-1.png differ
diff --git a/r-session-02_files/figure-html/unnamed-chunk-3-1.png b/r-session-02_files/figure-html/unnamed-chunk-3-1.png
index 723c5f3..d867486 100644
Binary files a/r-session-02_files/figure-html/unnamed-chunk-3-1.png and b/r-session-02_files/figure-html/unnamed-chunk-3-1.png differ
diff --git a/r-session-02_files/figure-html/unnamed-chunk-33-1.png b/r-session-02_files/figure-html/unnamed-chunk-33-1.png
index 1ebb2ef..b590e34 100644
Binary files a/r-session-02_files/figure-html/unnamed-chunk-33-1.png and b/r-session-02_files/figure-html/unnamed-chunk-33-1.png differ
diff --git a/r-session-02_files/figure-html/unnamed-chunk-35-1.png b/r-session-02_files/figure-html/unnamed-chunk-35-1.png
index 5501269..294e964 100644
Binary files a/r-session-02_files/figure-html/unnamed-chunk-35-1.png and b/r-session-02_files/figure-html/unnamed-chunk-35-1.png differ
diff --git a/r-session-02_files/figure-html/unnamed-chunk-62-1.png b/r-session-02_files/figure-html/unnamed-chunk-62-1.png
index c3b6cd1..8e617c7 100644
Binary files a/r-session-02_files/figure-html/unnamed-chunk-62-1.png and b/r-session-02_files/figure-html/unnamed-chunk-62-1.png differ
diff --git a/r-session-02_files/figure-html/unnamed-chunk-68-1.png b/r-session-02_files/figure-html/unnamed-chunk-68-1.png
index 8bfa465..af1da25 100644
Binary files a/r-session-02_files/figure-html/unnamed-chunk-68-1.png and b/r-session-02_files/figure-html/unnamed-chunk-68-1.png differ
diff --git a/r-session-02_files/figure-html/unnamed-chunk-71-1.png b/r-session-02_files/figure-html/unnamed-chunk-71-1.png
index 6cfc3e0..82267c0 100644
Binary files a/r-session-02_files/figure-html/unnamed-chunk-71-1.png and b/r-session-02_files/figure-html/unnamed-chunk-71-1.png differ
diff --git a/r-session-02_files/figure-html/unnamed-chunk-72-1.png b/r-session-02_files/figure-html/unnamed-chunk-72-1.png
index 8bdc722..6c093d9 100644
Binary files a/r-session-02_files/figure-html/unnamed-chunk-72-1.png and b/r-session-02_files/figure-html/unnamed-chunk-72-1.png differ
diff --git a/r-session-02_files/figure-html/unnamed-chunk-74-1.png b/r-session-02_files/figure-html/unnamed-chunk-74-1.png
index d37121d..59262e7 100644
Binary files a/r-session-02_files/figure-html/unnamed-chunk-74-1.png and b/r-session-02_files/figure-html/unnamed-chunk-74-1.png differ
diff --git a/r-session-02_files/figure-html/unnamed-chunk-76-1.png b/r-session-02_files/figure-html/unnamed-chunk-76-1.png
index 71f8f41..efb6be6 100644
Binary files a/r-session-02_files/figure-html/unnamed-chunk-76-1.png and b/r-session-02_files/figure-html/unnamed-chunk-76-1.png differ
diff --git a/r-session-02_files/figure-html/unnamed-chunk-78-1.png b/r-session-02_files/figure-html/unnamed-chunk-78-1.png
index 273f4ba..740a6e8 100644
Binary files a/r-session-02_files/figure-html/unnamed-chunk-78-1.png and b/r-session-02_files/figure-html/unnamed-chunk-78-1.png differ
diff --git a/r-session-02_files/figure-html/unnamed-chunk-8-1.png b/r-session-02_files/figure-html/unnamed-chunk-8-1.png
index 3383441..7c926fb 100644
Binary files a/r-session-02_files/figure-html/unnamed-chunk-8-1.png and b/r-session-02_files/figure-html/unnamed-chunk-8-1.png differ
diff --git a/r-session-02_files/figure-html/unnamed-chunk-80-1.png b/r-session-02_files/figure-html/unnamed-chunk-80-1.png
index 0d0dd58..4f5837a 100644
Binary files a/r-session-02_files/figure-html/unnamed-chunk-80-1.png and b/r-session-02_files/figure-html/unnamed-chunk-80-1.png differ
