From 0a9cddaa76cc451e0e78c8bb8344dfc63e3c63fa Mon Sep 17 00:00:00 2001 From: Irena Papst Date: Thu, 20 Mar 2025 16:28:31 -0400 Subject: [PATCH 1/3] rough first pass at age-structured SIR --- inst/starter_models/sir_age/README.Rmd | 176 +++++++++++++++++++++ inst/starter_models/sir_age/references.bib | 13 ++ inst/starter_models/sir_age/tmb.R | 35 ++++ 3 files changed, 224 insertions(+) create mode 100644 inst/starter_models/sir_age/README.Rmd create mode 100644 inst/starter_models/sir_age/references.bib create mode 100644 inst/starter_models/sir_age/tmb.R diff --git a/inst/starter_models/sir_age/README.Rmd b/inst/starter_models/sir_age/README.Rmd new file mode 100644 index 000000000..d28723d7f --- /dev/null +++ b/inst/starter_models/sir_age/README.Rmd @@ -0,0 +1,176 @@ +--- +title: "SIR with Age-based Transmission" +index_entry: "An SIR model with age-based transmission" +bibliography: ../../references.bib +link-citations: TRUE +author: Irena papst, Steve Walker +output: + github_document: + toc: true +editor_options: + chunk_output_type: console +--- + +```{r, include = FALSE} +knitr::opts_chunk$set( + collapse = TRUE, + warning = FALSE, + message = FALSE, + comment = "#>", + fig.path = "./figures/" +) +``` + +This is an extension of the [sir](https://github.com/canmod/macpan2/tree/main/inst/starter_models/sir) model that includes age-based transmission, based off of the model presented in [@Mistry2021]. + +# Packages Used and Settings + +The code in this article uses the following packages. + +```{r packages, message=FALSE, warning=FALSE} +library(ggplot2) +library(dplyr) +library(tidyr) +library(macpan2) +``` + +# Model Specification + +This model has been specified in the `sir_age` directory [here](https://github.com/canmod/macpan2/blob/main/inst/starter_models/sir_age/tmb.R) and is accessible from the `macpan2` model library (see [Example Models](https://canmod.github.io/macpan2/articles/example_models.html) for details). We can read in the model specification using the `mp_tmb_library` command. +```{r model_spec} +spec = mp_tmb_library( + "starter_models" + , "sir_age" + , package = "macpan2" +) +``` + +This specification can be used to draw the following flow diagram using code found in the [source for this article](https://github.com/canmod/macpan2/blob/main/inst/starter_models/sir_demog/README.Rmd). + +```{r diagram, echo = FALSE, eval = FALSE, fig.height = 2, fig.width = 8} +system.file("utilss", "box-drawing.R", package = "macpan2") |> source() +layout = mp_layout_paths(spec) +(layout + |> plot_flow_diagram(show_flow_names = TRUE) + |> draw_inflows(layout, show_labels = TRUE) + |> draw_outflows(layout, show_labels = TRUE, pattern_mutate = "^([a-z]*)(_.*)", pattern_replace = "\\1") +) +``` + + +# States + +| variable | description | +| -------- | --------------------------------- | +| S | Number of susceptible individuals | +| I | Number of infectious individuals | +| R | Number of recovered individuals | + +In this age structured model, each state is a vector with one entry per age group. + +The size of the population in each age groups is $N = S + I + R$. + +# Parameters + +| variable | description | +| -------- | ---------------------------- | +| $\tau$ | transmissibility (proportion of contacts with infecteds that yield infection) | +| $\gamma$ | per capita recovery rate | +| $M$ | matrix of contact rates between age groups, _e.g._ $M_{ij}$ is the rate of contact between susceptibles of age $i$ and infecteds of age $j$ | + +# Dynamics + +We assume new individuals (births) join the susceptible compartment, and individuals can leave the population (die) from any compartment. + +$$ +\begin{align*} +\frac{dS}{dt} &= - \lambda S \\ +\frac{dI}{dt} &= \lambda S - \gamma I \\ +\frac{dR}{dt} &= \gamma I +\end{align*} +$$ + +Here $\lambda = \tau M \frac{I}{N}$. + +# Simulation Example + +```{r plot_setup} +# variable labeller +matrix_labeller <- c( + "S" = "Susceptible", + "I" = "Infectious", + "R" = "Recovered" +) +row_labeller <- c( + "0" = "Children (ages 0-17)", + "1" = "Adults (ages 18-64)", + "2" = "Seniors (ages 65+)" +) + +plot_sim_age <- function(sim, matrix_labeller, row_labeller, title, subtitle){ + ggplot2::ggplot( + data = sim, + mapping = ggplot2::aes(x = time, y = value, colour = row) + ) + + ggplot2::geom_line() + + ggplot2::facet_wrap( + ~ matrix, + labeller = ggplot2::labeller( + matrix = matrix_labeller + ), + scales = "free_y") + + ggplot2::scale_colour_discrete( + labels = row_labeller + ) + + ggplot2::labs( + title = title, + subtitle = subtitle, + x = "Days from beginning of outbreak", + y = "Number of individuals" + ) + + ggplot2::theme( + legend.position = "bottom" + , legend.title = ggplot2::element_blank() + ) +} +``` + +```{r plot_sim} +# set parameters +N <- rep(1e5, 3) # total population by age +I0 <- 100 # initial infecteds + +M <- matrix(c(10, 1, 1, 1, 5, 1, 1, 1, 2), nrow = 3) # contact matrix in units of avg number of contacts/day (symmetric contact rates because we will pick two age groups of the same size) +gamma <- 1/5 # 5 days avg infectious period; same gamma across age groups + +R0 <- 2 +# spectral radius of contact matrix for transmissibility calculation +rho <- max(abs(eigen(M)$values)) +# calculate transmissibility given R0, gamma, M +tau <- R0*gamma/rho # transmissibility (proportion of contacts with infectious individuals that yield infection) + +# create simulator +simulator <- mp_simulator( + model = spec, + time_steps = 150 + , outputs = mp_state_vars(spec) +) + +# generate trajectory +sim <- mp_trajectory(simulator) |> + dplyr::mutate( + matrix = factor(matrix, levels = mp_state_vars(spec)) + , row = forcats::as_factor(row)) # natural order of state variables + +# visualize trajectory +plot_sim_age( + sim, + matrix_labeller, + row_labeller, + title = "Simulation of the age-structured SIR model with two age groups", + + subtitle = paste("Infection is seeded with", I0, "infectious children in a total population of", scales::label_number(scale_cut = scales::cut_short_scale())(sum(N))) +) +``` + +# References diff --git a/inst/starter_models/sir_age/references.bib b/inst/starter_models/sir_age/references.bib new file mode 100644 index 000000000..ce5f6f045 --- /dev/null +++ b/inst/starter_models/sir_age/references.bib @@ -0,0 +1,13 @@ +@article{Mistry2021, + title = {Inferring high-resolution human mixing patterns for disease modeling}, + volume = {12}, + ISSN = {2041-1723}, + url = {http://dx.doi.org/10.1038/s41467-020-20544-y}, + DOI = {10.1038/s41467-020-20544-y}, + number = {1}, + journal = {Nature Communications}, + publisher = {Springer Science and Business Media LLC}, + author = {Mistry, Dina and Litvinova, Maria and Pastore y Piontti, Ana and Chinazzi, Matteo and Fumanelli, Laura and Gomes, Marcelo F. C. and Haque, Syed A. and Liu, Quan-Hui and Mu, Kunpeng and Xiong, Xinyue and Halloran, M. Elizabeth and Longini, Ira M. and Merler, Stefano and Ajelli, Marco and Vespignani, Alessandro}, + year = {2021}, + month = jan +} \ No newline at end of file diff --git a/inst/starter_models/sir_age/tmb.R b/inst/starter_models/sir_age/tmb.R new file mode 100644 index 000000000..280186fc1 --- /dev/null +++ b/inst/starter_models/sir_age/tmb.R @@ -0,0 +1,35 @@ +library(macpan2) + +# SIR with age, where age structure only affects transmission, as +# defined here: https://www.nature.com/articles/s41467-020-20544-y#Sec7 + +# define the model +spec = mp_tmb_model_spec( + before = list( + S ~ N - I - R # calculate population size before the simulation loop begins to avoid having to specify a value for it by hand in the defaults list + ), + during = list( + # force of infection + lambda ~ tau * M %*% (I / N), + # infection + mp_per_capita_flow( + from = "S", to = "I", + rate = "lambda", + abs_rate = "infection" + ), + # recovery + mp_per_capita_flow( + from = "I", to = "R", + rate = "gamma", + abs_rate = "recovery" + ) + ), + default = list( + tau = tau, + M = M, + gamma = gamma, # 5 days avg infectious period; same gamma across age groups + N = N, + I = c(I0, 0, 0), # seed infection in age group 1 + R = rep(0, 3) + ) +) From c3f58e4f444c4385295907c7de316c394cd9f57f Mon Sep 17 00:00:00 2001 From: Irena Papst Date: Thu, 20 Mar 2025 16:30:41 -0400 Subject: [PATCH 2/3] Update README.Rmd --- inst/starter_models/sir_age/README.Rmd | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/inst/starter_models/sir_age/README.Rmd b/inst/starter_models/sir_age/README.Rmd index d28723d7f..22f022652 100644 --- a/inst/starter_models/sir_age/README.Rmd +++ b/inst/starter_models/sir_age/README.Rmd @@ -3,7 +3,7 @@ title: "SIR with Age-based Transmission" index_entry: "An SIR model with age-based transmission" bibliography: ../../references.bib link-citations: TRUE -author: Irena papst, Steve Walker +author: Irena Papst, Steve Walker output: github_document: toc: true From 5fc3cb48f07a24d91de4770f2cb3e88ed99de085 Mon Sep 17 00:00:00 2001 From: Irena Papst Date: Fri, 21 Mar 2025 07:44:39 -0400 Subject: [PATCH 3/3] tweaks to README --- inst/starter_models/sir_age/README.Rmd | 120 ++++++++++++------------- 1 file changed, 56 insertions(+), 64 deletions(-) diff --git a/inst/starter_models/sir_age/README.Rmd b/inst/starter_models/sir_age/README.Rmd index 22f022652..209cba5a2 100644 --- a/inst/starter_models/sir_age/README.Rmd +++ b/inst/starter_models/sir_age/README.Rmd @@ -1,6 +1,6 @@ --- -title: "SIR with Age-based Transmission" -index_entry: "An SIR model with age-based transmission" +title: "Age-stratified SIR" +index_entry: "An age-stratified SIR model" bibliography: ../../references.bib link-citations: TRUE author: Irena Papst, Steve Walker @@ -21,7 +21,7 @@ knitr::opts_chunk$set( ) ``` -This is an extension of the [sir](https://github.com/canmod/macpan2/tree/main/inst/starter_models/sir) model that includes age-based transmission, based off of the model presented in [@Mistry2021]. +This is an expansion of the [sir](https://github.com/canmod/macpan2/tree/main/inst/starter_models/sir) model to include age-stratification; each state variable is subdivided by age, and transmission is governed by a age-based contact matrix, as described in [@Mistry2021]. # Packages Used and Settings @@ -30,7 +30,8 @@ The code in this article uses the following packages. ```{r packages, message=FALSE, warning=FALSE} library(ggplot2) library(dplyr) -library(tidyr) +library(forcats) +library(scales) library(macpan2) ``` @@ -45,7 +46,7 @@ spec = mp_tmb_library( ) ``` -This specification can be used to draw the following flow diagram using code found in the [source for this article](https://github.com/canmod/macpan2/blob/main/inst/starter_models/sir_demog/README.Rmd). +This specification can be used to draw the following flow diagram using code found in the [source for this article](https://github.com/canmod/macpan2/blob/main/inst/starter_models/sir_age/README.Rmd). ```{r diagram, echo = FALSE, eval = FALSE, fig.height = 2, fig.width = 8} system.file("utilss", "box-drawing.R", package = "macpan2") |> source() @@ -66,21 +67,21 @@ layout = mp_layout_paths(spec) | I | Number of infectious individuals | | R | Number of recovered individuals | -In this age structured model, each state is a vector with one entry per age group. +Each state variable is a vector with one entry per age group. -The size of the population in each age groups is $N = S + I + R$. +The size of the population in each age groups is given by the components of $N = S + I + R$. # Parameters | variable | description | | -------- | ---------------------------- | -| $\tau$ | transmissibility (proportion of contacts with infecteds that yield infection) | +| $\tau$ | transmissibility, _i.e._, the proportion of contacts with infecteds that yield infection | | $\gamma$ | per capita recovery rate | | $M$ | matrix of contact rates between age groups, _e.g._ $M_{ij}$ is the rate of contact between susceptibles of age $i$ and infecteds of age $j$ | -# Dynamics +$\tau$ and $\gamma$ can be specified by age (as vectors), though in the example below we use a scalar value that applies to all age groups. -We assume new individuals (births) join the susceptible compartment, and individuals can leave the population (die) from any compartment. +# Dynamics $$ \begin{align*} @@ -94,8 +95,25 @@ Here $\lambda = \tau M \frac{I}{N}$. # Simulation Example -```{r plot_setup} -# variable labeller +```{r plot_sim} +# set numerical values +N <- rep(1e5, 3) # total population by age +I0 <- 100 # initial infecteds + +M <- matrix(c( + 10, 1, 1, + 1, 5, 1, + 1, 1, 2 + ), nrow = 3) # contact matrix in units of avg number of contacts/person/day (symmetric contact rates because we will pick two age groups of the same size) +gamma <- 1/5 # 5 days avg infectious period; same gamma across age groups + +R0 <- 2 +# spectral radius of contact matrix for transmissibility calculation +rho <- max(abs(eigen(M)$values)) +# calculate transmissibility given R0, gamma, M +tau <- R0*gamma/rho # transmissibility (proportion of contacts with infectious individuals that yield infection) + +# variable labellers matrix_labeller <- c( "S" = "Susceptible", "I" = "Infectious", @@ -107,69 +125,43 @@ row_labeller <- c( "2" = "Seniors (ages 65+)" ) -plot_sim_age <- function(sim, matrix_labeller, row_labeller, title, subtitle){ - ggplot2::ggplot( - data = sim, +# simulate & visualize +(spec + # simulate + |> mp_simulator( + time_steps = 150 + , outputs = mp_state_vars(spec) + ) + |> mp_trajectory() + # enforce natural order of variables for plotting + |> dplyr::mutate( + matrix = factor(matrix, levels = mp_state_vars(spec)) + , row = forcats::as_factor(row) + ) + # visualize trajectory + |> ggplot2::ggplot( mapping = ggplot2::aes(x = time, y = value, colour = row) - ) + - ggplot2::geom_line() + - ggplot2::facet_wrap( + ) + + ggplot2::geom_line() + + + ggplot2::facet_wrap( ~ matrix, labeller = ggplot2::labeller( matrix = matrix_labeller ), - scales = "free_y") + - ggplot2::scale_colour_discrete( + scales = "free_y") + + ggplot2::scale_colour_discrete( labels = row_labeller - ) + - ggplot2::labs( - title = title, - subtitle = subtitle, + ) + + ggplot2::labs( + title = "Simulation of the age-structured SIR model with two age groups", + subtitle = paste("Infection is seeded with", I0, "infectious children in a total population of", scales::label_number(scale_cut = scales::cut_short_scale())(sum(N))), x = "Days from beginning of outbreak", y = "Number of individuals" - ) + - ggplot2::theme( + ) + + ggplot2::theme( legend.position = "bottom" , legend.title = ggplot2::element_blank() ) -} -``` - -```{r plot_sim} -# set parameters -N <- rep(1e5, 3) # total population by age -I0 <- 100 # initial infecteds - -M <- matrix(c(10, 1, 1, 1, 5, 1, 1, 1, 2), nrow = 3) # contact matrix in units of avg number of contacts/day (symmetric contact rates because we will pick two age groups of the same size) -gamma <- 1/5 # 5 days avg infectious period; same gamma across age groups - -R0 <- 2 -# spectral radius of contact matrix for transmissibility calculation -rho <- max(abs(eigen(M)$values)) -# calculate transmissibility given R0, gamma, M -tau <- R0*gamma/rho # transmissibility (proportion of contacts with infectious individuals that yield infection) - -# create simulator -simulator <- mp_simulator( - model = spec, - time_steps = 150 - , outputs = mp_state_vars(spec) -) - -# generate trajectory -sim <- mp_trajectory(simulator) |> - dplyr::mutate( - matrix = factor(matrix, levels = mp_state_vars(spec)) - , row = forcats::as_factor(row)) # natural order of state variables - -# visualize trajectory -plot_sim_age( - sim, - matrix_labeller, - row_labeller, - title = "Simulation of the age-structured SIR model with two age groups", - - subtitle = paste("Infection is seeded with", I0, "infectious children in a total population of", scales::label_number(scale_cut = scales::cut_short_scale())(sum(N))) ) ```