From ad4b5daeeb561c8bde54b880eefcb07f882f14ca Mon Sep 17 00:00:00 2001 From: Callum Arnold <24391445+arnold-c@users.noreply.github.com> Date: Mon, 7 Jul 2025 20:02:23 -0400 Subject: [PATCH] Update README with description of the project --- README.md | 118 +++++++++++++++++- _book/r-session-04.html | 2 +- _book/search.json | 2 +- .../r-session-04/execute-results/html.json | 4 +- 4 files changed, 117 insertions(+), 9 deletions(-) diff --git a/README.md b/README.md index 7d7059b..5408440 100644 --- a/README.md +++ b/README.md @@ -1,23 +1,131 @@ # SISMID Module 2 Materials (2025) + +**Mathematical Models of Infectious Diseases** + +This repository contains the educational materials for the [2025 Summer Institute in Statistics and Modeling in Infectious Diseases (SISMID) Module 2: Mathematical Models of Infectious Diseases](https://sph.emory.edu/SISMID/modules/math-models/index.html). The materials are presented as an interactive Quarto book website with lecture notes, hands-on R exercises, and comprehensive tutorials for epidemiological modeling. + ## About This Project +This educational resource provides a comprehensive introduction to mathematical modeling of infectious diseases, covering fundamental concepts from basic SIR dynamics to advanced topics like stochastic models and parameter estimation. +The materials are designed for a 3-day intensive workshop format but can also be used for self-paced learning. + +### Content Overview + +The materials are organized into three main sections: + +**Pre-Requisites**: Setup guides and foundational knowledge +- R and RStudio installation +- Essential R programming concepts +- Project organization best practices + +**Day 1**: Introduction and Basic Models +- Introduction to modeling concepts +- Basic SIR dynamics +- Herd immunity and vaccination strategies + +**Day 2**: Advanced Modeling Concepts +- Heterogeneity in models +- Age-structured models +- Parameter estimation techniques + +**Day 3**: Stochastic Models and Applications +- Stochastic modeling approaches +- Advanced R exercises and applications + ### Who Is This For & Pre-Requisite Knowledge? +This material is designed for: +- Graduate students and researchers in epidemiology, public health, or related fields +- Practitioners working in infectious disease modeling +- Anyone interested in learning mathematical approaches to understanding disease transmission + +**Prerequisites:** +- Basic understanding of calculus and differential equations +- Familiarity with R programming (see the ["Just Enough R" section](https://sismid2025.callumarnold.com/just-enough-r) if you're new to R) +- Working installation of R and RStudio +- Approximately 30-60 minutes for setup and prerequisite review + +## Getting Started + +### Download and Build the Website + +1. **Clone the repository:** + ```bash + git clone https://github.com/arnold-c/SISMID-Module-02_2025.git + cd SISMID-Module-02_2025 + ``` + +2. **Install R dependencies:** + This project uses `renv` for reproducible package management. In R: + ```r + # Install renv if you don't have it + install.packages("renv") + + # Restore the project library + renv::restore() + ``` + +3. **Install Quarto:** + Download and install Quarto from [https://quarto.org/docs/get-started/](https://quarto.org/docs/get-started/) + +4. **Build the website:** + ```bash + # Render the entire book + quarto render + + # Or preview with live reload during development + quarto preview + ``` + +5. **View the website:** + Open `_book/index.html` in your browser, or use the local server from `quarto preview` + +### Development Workflow + +- **Live preview**: `quarto preview` - starts a local server with automatic reload +- **Render single file**: `quarto render filename.qmd` +- **Check project**: `quarto check` +- **Update R packages**: `renv::snapshot()` after adding new packages ## Built With -- [Quarto](https://quarto.org) +- [Quarto](https://quarto.org) - Scientific and technical publishing system +- [R](https://www.r-project.org/) - Statistical computing and graphics +- [renv](https://rstudio.github.io/renv/) - Reproducible R package management +- Key R packages: `{tidyverse}`, `{deSolve`}, `{ggplot2`}, `{gt`}, and more (see `renv.lock`) ## Contact -You can contact me via my email: "callum *at* callumarnold *dot* com". +You can contact me via my email: "contact *at* callumarnold *dot* com". ## Contributing -If you see any issues, please open an Issue or Discussion on the book's GitHub page. -Or, if you know Git, you can open a Pull Request to fix the issue directly, which is even better! +We welcome contributions to improve these educational materials! Here are several ways you can help: + +### Reporting Issues +- **Found a bug or error?** Open an [Issue](https://github.com/arnold-c/SISMID-Module-02_2025/issues) describing the problem +- **Have a suggestion?** Start a [Discussion](https://github.com/arnold-c/SISMID-Module-02_2025/discussions) to share your ideas +- **Unclear instructions?** Let us know what could be explained better (open a [Discussion](https://github.com/arnold-c/SISMID-Module-02_2025/discussions)) + +### Contributing Code or Content +1. **Fork the repository** and create a feature branch +2. **Make your changes** following the existing style and structure + - The R code uses the [Air](https://posit-dev.github.io/air/) formatter for the R code chunks. Please follow the instructions on their website for the best way to employ the formatter in your development environment. +3. **Test your changes** by rendering the affected pages with `quarto render` +4. **Submit a Pull Request** with a clear description of your changes + +### Content Guidelines +- Follow the existing Quarto markdown format and callout patterns +- Test all R code chunks to ensure they execute correctly +- Maintain consistency with the educational tone and level +- Include appropriate citations for new references + +### Development Setup for Contributors +- Ensure you can successfully run `quarto preview` and `renv::restore()` +- Check that your changes don't break existing functionality +- Follow the project's code style and organization patterns -If you are interested in contributing text, please email me and we can make it happen! +**Interested in contributing substantial content?** Please reach out to me first by opening a discussion on this GitHub page (and @ me) so we can coordinate efforts and ensure your contributions align with the overall educational goals. ## Acknowledgements diff --git a/_book/r-session-04.html b/_book/r-session-04.html index 0455339..66df953 100644 --- a/_book/r-session-04.html +++ b/_book/r-session-04.html @@ -2179,7 +2179,7 @@ R <- 1 # number recovered time <- 0 -beta <- .5 # transmission rate +beta <- 0.5 # transmission rate gamma <- 1 / 7 # recovery rate ``` diff --git a/_book/search.json b/_book/search.json index bdbfadd..7978bc1 100644 --- a/_book/search.json +++ b/_book/search.json @@ -698,7 +698,7 @@ "href": "r-session-04.html#the-gillespie-algorithm", "title": "\n11  R Session 04\n", "section": "\n11.3 The Gillespie Algorithm", - "text": "11.3 The Gillespie Algorithm\nThe Gillespie algorithm is the most explicit translation of the ODE form of the SIR model into a stochastic model. It achieves this by noting that ODE-based models are written in terms of the rates at which events occur. At any given point in time, the rate of all events in the ODE is known; what known, is which of the possible events will happen first. Note that even though we expect that the next thing to happen will be the thing that happens with the highest rate, it’s possible that another thing may happen first by random chance. And then, since the rates in the SIR model are dependent on the value of the states (e.g. new infections depend on \\(\\beta\\) and S and I) any change in the states then changes the rates and the likelihood of what will happen next.\nThe Gillespie algorithm proceeds by 1) calculating all the current rates, 2) randomly drawing exponential random variables that equate to the time until each event happens (recall we talked about the relationship between rate and time in the lectures), then 3) comparing those times and assuming that the next event to occur is the one that had the smallest randomly drawn time. Then you increment the states; e.g. an infection increases I by 1, reduces S by 1, and doesn’t change R. Importantly, you then increment time forward by a step equal to the time until first event occurred. Then you recalculate the rates and randomly draw times, etc, and keep doing this over and over again. For this SIR model, that keeps happening until you run out of infected individuals and there are no new events that can happen.\n\nCode####################################################################\n# Parameters and initial conditions #\n####################################################################\n\nS <- 998 # number susceptible\nI <- 1 # number infected\nR <- 1 # number recovered\ntime <- 0\n\nbeta <- .5 # transmission rate\ngamma <- 1 / 7 # recovery rate\n\n\n\nCode####################################################################\n# Gillespie Step #\n####################################################################\ngil_step <- function(SIR, beta, gamma) {\n # SIR is a vector containing 4 elements\n # S = scalar number of susceptibles now\n # I = scalar number of infecteds now\n # R = scalar number of recovereds now\n # time = current time\n # beta = transmission rate\n # gamma = recovery rate\n\n # draw two random exponential variables\n times <- rexp(\n 2,\n c(\n # the first is the rate of new infections: S*Beta*I/N\n beta * SIR[1] * SIR[2] / sum(SIR[1:3]),\n # the second is the rate of new recoveries: I*gamma\n SIR[2] * gamma\n )\n )\n\n return(list(\n # which.min identifies which of the two random variables is smallest,\n # and thus the first to happen: this is the transition that we'll make\n change_state = which.min(times),\n # this identifies how much time elapsed before the transition occurred,\n # so we can increment time forward\n time_step = min(times)\n ))\n}\n\n\nThe gil_step() function increments the states forward by 1 event. Which event happens first, and the time it takes for that event to happen are both random variables. We then need to do this many times.\n\nCode####################################################################\n# Simulate over time #\n####################################################################\n# here we set the random seed. This isn't necessary in general, but it allows us to write a \"random\" simulation\nset.seed(101)\ncounter <- 0 # set counter at 0\nwhile (all(I > 0)) {\n # continue until I is depleted\n counter <- counter + 1\n # counter for number of transitions: we don't know how many transitions\n # will happen until the simulation is over\n #\n # current SIR states\n sir_tmp <- c(S[counter], I[counter], R[counter], time[counter])\n step <- gil_step(sir_tmp, beta, gamma)\n if (step$change_state == 1) {\n # if transition is an infection, reduce S and increase I\n sir_tmp[1] <- sir_tmp[1] - 1\n sir_tmp[2] <- sir_tmp[2] + 1\n }\n if (step$change_state == 2) {\n # if transition is an recovery, reduce I and increase R\n sir_tmp[2] <- sir_tmp[2] - 1\n sir_tmp[3] <- sir_tmp[3] + 1\n }\n sir_tmp[4] <- sir_tmp[4] + step$time_step # increment time\n\n # Append changes\n S <- c(S, sir_tmp[1])\n I <- c(I, sir_tmp[2])\n R <- c(R, sir_tmp[3])\n time <- c(time, sir_tmp[4])\n\n # cat(S[counter],\\\"-\\\",I[counter],\\\"-\\\",R[counter],\\\"-\\\",time[counter], \\\".\\\\n\\\")\n # this prints the output as it goes\n # reset the seed so that every subsequent simulation IS random\n set.seed(NULL)\n}\n\n\nWe can now plot the one realization of this stochastic outbreak. Notice that it has the same general shape as the deterministic outbreak, but is no longer smooth because each individual event over time happens randomly.\n\nCode# Reuse the plotting code from above\n\n# Turn the output from gil_step() into a tibble (dataframe)\n# so we can manipulate and plot it easily\n# put elements in a data frame\ngil_df <- tibble(time, S, I, R) %>%\n # Convert all columns to numeric (they are currently type\n # deSolve so will produce warnings when plotting etc)\n mutate(\n # Rather than repeatedly type the same function for every\n # column, use the across() function to apply the function\n # to a selection of columns\n across(\n # The cols argument takes a selection of columns to apply\n # a function to. Here, we want to apply the as.numeric()\n # function to all columns, so we use the function\n # everything() to select all columns.\n .cols = everything(),\n .fns = as.numeric\n )\n ) %>%\n # Convert the dataframe from wide to long format, so we have a\n # column for the time, a column for the state, and a column\n # for the proportion of the population in that state at that\n # time\n pivot_longer(\n # Don't pivot the time column\n cols = -c(time),\n names_to = \"state\",\n values_to = \"number\"\n ) %>%\n # Update the state column to be a factor, so the plot will\n # show the states in the correct order\n mutate(state = factor(state, levels = c(\"S\", \"I\", \"R\")))\n\nsir_colors <- c(S = \"#1f77b4\", I = \"#ff7f0e\", R = \"#FF3851\")\n\nggplot(gil_df, aes(x = time, y = number, color = state)) +\n geom_line(linewidth = 1.5) +\n scale_color_manual(values = sir_colors) +\n labs(\n x = \"Time\",\n y = \"Number\",\n color = \"State\"\n ) +\n theme(legend.position = \"top\")\n\n\n\n\n\n\n\nBecause each realization is stochastic, we need to generate many runs to see the general behavior. The code below runs 10 iterations. Generally, this would be considered a very small number of iterations. But even for this very small model, running 100 or more means you’ll be waiting for output. Note that this code is designed to be transparent not to be fast. Making these run fast is beyond the scope of this assignment.\n\nCode####################################################################\n# Simulate over time #\nnum_iterations <- 10 # number of realizations to simulate\n\nfor (iter in 1:num_iterations) {\n ####################################################################\n # Parameters and initial conditions #\n ####################################################################\n\n S <- 998 # number susceptible\n I <- 1 # number infected\n R <- 1 # number recovered\n time <- 0\n # initialize iteration counter\n iteration <- iter\n\n beta <- .5 # transmission rate\n gamma <- 1 / 7 # recovery rate\n\n ####################################################################\n # Simulate over time #\n ####################################################################\n counter <- 0 # set counter at 0\n while (all(I > 0)) {\n # continue until I is depleted\n counter <- counter + 1\n # counter for number of transitions: we don't know how many transitions\n # will happen until the simulation is over\n\n # current SIR states\n sir_tmp <- c(S[counter], I[counter], R[counter], time[counter])\n step <- gil_step(sir_tmp, beta, gamma)\n if (step$change_state == 1) {\n # if transition is an infection, reduce S and increase I\n sir_tmp[1] <- sir_tmp[1] - 1\n sir_tmp[2] <- sir_tmp[2] + 1\n }\n if (step$change_state == 2) {\n # if transition is an recovery, reduce I and increase R\n sir_tmp[2] <- sir_tmp[2] - 1\n sir_tmp[3] <- sir_tmp[3] + 1\n }\n # increment time\n sir_tmp[4] <- sir_tmp[4] + step$time_step\n\n # Append changes\n S <- c(S, sir_tmp[1])\n I <- c(I, sir_tmp[2])\n R <- c(R, sir_tmp[3])\n time <- c(time, sir_tmp[4])\n iteration <- c(iteration, iter)\n\n # cat(S[counter],\\\"-\\\",I[counter],\\\"-\\\",R[counter],\\\"-\\\",time[counter], \\\".\\\\n\\\")\n # this prints the output as it goes\n }\n if (iter == 1) {\n sir_gil_storage <- tibble(S, I, R, time, iteration)\n }\n if (iter > 1) {\n sir_gil_storage <- bind_rows(\n sir_gil_storage,\n tibble(S, I, R, time, iteration)\n )\n }\n}\n\n\n\nCode# Reuse the plotting code from above\n\n# Turn the output from gil_step() into a tibble (dataframe)\n# so we can manipulate and plot it easily\n# put elements in a data frame\ngil_df <- sir_gil_storage %>%\n # Convert all columns to numeric (they are currently type\n # deSolve so will produce warnings when plotting etc)\n mutate(\n # Rather than repeatedly type the same function for every\n # column, use the across() function to apply the function\n # to a selection of columns\n across(\n # The cols argument takes a selection of columns to apply\n # a function to. Here, we want to apply the as.numeric()\n # function to all columns, so we use the function\n # everything() to select all columns.\n .cols = everything(),\n .fns = as.numeric\n )\n ) %>%\n # Convert the dataframe from wide to long format, so we have a\n # column for the time, a column for the state, and a column\n # for the proportion of the population in that state at that\n # time\n pivot_longer(\n # Don't pivot the time column\n cols = -c(time, iteration),\n names_to = \"state\",\n values_to = \"number\"\n ) %>%\n # Update the state column to be a factor, so the plot will\n # show the states in the correct order\n mutate(state = factor(state, levels = c(\"S\", \"I\", \"R\")))\n\nsir_colors <- c(S = \"#1f77b4\", I = \"#ff7f0e\", R = \"#FF3851\")\n\nggplot(\n gil_df,\n aes(\n x = time,\n y = number,\n color = state,\n group = interaction(iteration, state)\n )\n) +\n geom_line(linewidth = 1.5, alpha = .1) +\n scale_color_manual(values = sir_colors) +\n labs(\n x = \"Time\",\n y = \"Number\",\n color = \"State\"\n ) +\n guides(color = guide_legend(override.aes = list(alpha = 1))) +\n theme(legend.position = \"top\")\n\n\n\n\n\n\n\nThe Gillespie algorithm doesn’t cut any corners relative to the ODE model, but that comes at the cost of computational efficiency. For every event that happens (e.g. an infection) you also have to generate a random draw (e.g. a recovery) that you don’t use, except for comparison. And the bigger your population, the more possible events can happen. So the bigger the model (e.g. adding exposed or vaccinated classes, or heterogeneity, add transitions) and the bigger the population, the more calculation you need and and therefore the slower the model.", + "text": "11.3 The Gillespie Algorithm\nThe Gillespie algorithm is the most explicit translation of the ODE form of the SIR model into a stochastic model. It achieves this by noting that ODE-based models are written in terms of the rates at which events occur. At any given point in time, the rate of all events in the ODE is known; what known, is which of the possible events will happen first. Note that even though we expect that the next thing to happen will be the thing that happens with the highest rate, it’s possible that another thing may happen first by random chance. And then, since the rates in the SIR model are dependent on the value of the states (e.g. new infections depend on \\(\\beta\\) and S and I) any change in the states then changes the rates and the likelihood of what will happen next.\nThe Gillespie algorithm proceeds by 1) calculating all the current rates, 2) randomly drawing exponential random variables that equate to the time until each event happens (recall we talked about the relationship between rate and time in the lectures), then 3) comparing those times and assuming that the next event to occur is the one that had the smallest randomly drawn time. Then you increment the states; e.g. an infection increases I by 1, reduces S by 1, and doesn’t change R. Importantly, you then increment time forward by a step equal to the time until first event occurred. Then you recalculate the rates and randomly draw times, etc, and keep doing this over and over again. For this SIR model, that keeps happening until you run out of infected individuals and there are no new events that can happen.\n\nCode####################################################################\n# Parameters and initial conditions #\n####################################################################\n\nS <- 998 # number susceptible\nI <- 1 # number infected\nR <- 1 # number recovered\ntime <- 0\n\nbeta <- 0.5 # transmission rate\ngamma <- 1 / 7 # recovery rate\n\n\n\nCode####################################################################\n# Gillespie Step #\n####################################################################\ngil_step <- function(SIR, beta, gamma) {\n # SIR is a vector containing 4 elements\n # S = scalar number of susceptibles now\n # I = scalar number of infecteds now\n # R = scalar number of recovereds now\n # time = current time\n # beta = transmission rate\n # gamma = recovery rate\n\n # draw two random exponential variables\n times <- rexp(\n 2,\n c(\n # the first is the rate of new infections: S*Beta*I/N\n beta * SIR[1] * SIR[2] / sum(SIR[1:3]),\n # the second is the rate of new recoveries: I*gamma\n SIR[2] * gamma\n )\n )\n\n return(list(\n # which.min identifies which of the two random variables is smallest,\n # and thus the first to happen: this is the transition that we'll make\n change_state = which.min(times),\n # this identifies how much time elapsed before the transition occurred,\n # so we can increment time forward\n time_step = min(times)\n ))\n}\n\n\nThe gil_step() function increments the states forward by 1 event. Which event happens first, and the time it takes for that event to happen are both random variables. We then need to do this many times.\n\nCode####################################################################\n# Simulate over time #\n####################################################################\n# here we set the random seed. This isn't necessary in general, but it allows us to write a \"random\" simulation\nset.seed(101)\ncounter <- 0 # set counter at 0\nwhile (all(I > 0)) {\n # continue until I is depleted\n counter <- counter + 1\n # counter for number of transitions: we don't know how many transitions\n # will happen until the simulation is over\n #\n # current SIR states\n sir_tmp <- c(S[counter], I[counter], R[counter], time[counter])\n step <- gil_step(sir_tmp, beta, gamma)\n if (step$change_state == 1) {\n # if transition is an infection, reduce S and increase I\n sir_tmp[1] <- sir_tmp[1] - 1\n sir_tmp[2] <- sir_tmp[2] + 1\n }\n if (step$change_state == 2) {\n # if transition is an recovery, reduce I and increase R\n sir_tmp[2] <- sir_tmp[2] - 1\n sir_tmp[3] <- sir_tmp[3] + 1\n }\n sir_tmp[4] <- sir_tmp[4] + step$time_step # increment time\n\n # Append changes\n S <- c(S, sir_tmp[1])\n I <- c(I, sir_tmp[2])\n R <- c(R, sir_tmp[3])\n time <- c(time, sir_tmp[4])\n\n # cat(S[counter],\\\"-\\\",I[counter],\\\"-\\\",R[counter],\\\"-\\\",time[counter], \\\".\\\\n\\\")\n # this prints the output as it goes\n # reset the seed so that every subsequent simulation IS random\n set.seed(NULL)\n}\n\n\nWe can now plot the one realization of this stochastic outbreak. Notice that it has the same general shape as the deterministic outbreak, but is no longer smooth because each individual event over time happens randomly.\n\nCode# Reuse the plotting code from above\n\n# Turn the output from gil_step() into a tibble (dataframe)\n# so we can manipulate and plot it easily\n# put elements in a data frame\ngil_df <- tibble(time, S, I, R) %>%\n # Convert all columns to numeric (they are currently type\n # deSolve so will produce warnings when plotting etc)\n mutate(\n # Rather than repeatedly type the same function for every\n # column, use the across() function to apply the function\n # to a selection of columns\n across(\n # The cols argument takes a selection of columns to apply\n # a function to. Here, we want to apply the as.numeric()\n # function to all columns, so we use the function\n # everything() to select all columns.\n .cols = everything(),\n .fns = as.numeric\n )\n ) %>%\n # Convert the dataframe from wide to long format, so we have a\n # column for the time, a column for the state, and a column\n # for the proportion of the population in that state at that\n # time\n pivot_longer(\n # Don't pivot the time column\n cols = -c(time),\n names_to = \"state\",\n values_to = \"number\"\n ) %>%\n # Update the state column to be a factor, so the plot will\n # show the states in the correct order\n mutate(state = factor(state, levels = c(\"S\", \"I\", \"R\")))\n\nsir_colors <- c(S = \"#1f77b4\", I = \"#ff7f0e\", R = \"#FF3851\")\n\nggplot(gil_df, aes(x = time, y = number, color = state)) +\n geom_line(linewidth = 1.5) +\n scale_color_manual(values = sir_colors) +\n labs(\n x = \"Time\",\n y = \"Number\",\n color = \"State\"\n ) +\n theme(legend.position = \"top\")\n\n\n\n\n\n\n\nBecause each realization is stochastic, we need to generate many runs to see the general behavior. The code below runs 10 iterations. Generally, this would be considered a very small number of iterations. But even for this very small model, running 100 or more means you’ll be waiting for output. Note that this code is designed to be transparent not to be fast. Making these run fast is beyond the scope of this assignment.\n\nCode####################################################################\n# Simulate over time #\nnum_iterations <- 10 # number of realizations to simulate\n\nfor (iter in 1:num_iterations) {\n ####################################################################\n # Parameters and initial conditions #\n ####################################################################\n\n S <- 998 # number susceptible\n I <- 1 # number infected\n R <- 1 # number recovered\n time <- 0\n # initialize iteration counter\n iteration <- iter\n\n beta <- .5 # transmission rate\n gamma <- 1 / 7 # recovery rate\n\n ####################################################################\n # Simulate over time #\n ####################################################################\n counter <- 0 # set counter at 0\n while (all(I > 0)) {\n # continue until I is depleted\n counter <- counter + 1\n # counter for number of transitions: we don't know how many transitions\n # will happen until the simulation is over\n\n # current SIR states\n sir_tmp <- c(S[counter], I[counter], R[counter], time[counter])\n step <- gil_step(sir_tmp, beta, gamma)\n if (step$change_state == 1) {\n # if transition is an infection, reduce S and increase I\n sir_tmp[1] <- sir_tmp[1] - 1\n sir_tmp[2] <- sir_tmp[2] + 1\n }\n if (step$change_state == 2) {\n # if transition is an recovery, reduce I and increase R\n sir_tmp[2] <- sir_tmp[2] - 1\n sir_tmp[3] <- sir_tmp[3] + 1\n }\n # increment time\n sir_tmp[4] <- sir_tmp[4] + step$time_step\n\n # Append changes\n S <- c(S, sir_tmp[1])\n I <- c(I, sir_tmp[2])\n R <- c(R, sir_tmp[3])\n time <- c(time, sir_tmp[4])\n iteration <- c(iteration, iter)\n\n # cat(S[counter],\\\"-\\\",I[counter],\\\"-\\\",R[counter],\\\"-\\\",time[counter], \\\".\\\\n\\\")\n # this prints the output as it goes\n }\n if (iter == 1) {\n sir_gil_storage <- tibble(S, I, R, time, iteration)\n }\n if (iter > 1) {\n sir_gil_storage <- bind_rows(\n sir_gil_storage,\n tibble(S, I, R, time, iteration)\n )\n }\n}\n\n\n\nCode# Reuse the plotting code from above\n\n# Turn the output from gil_step() into a tibble (dataframe)\n# so we can manipulate and plot it easily\n# put elements in a data frame\ngil_df <- sir_gil_storage %>%\n # Convert all columns to numeric (they are currently type\n # deSolve so will produce warnings when plotting etc)\n mutate(\n # Rather than repeatedly type the same function for every\n # column, use the across() function to apply the function\n # to a selection of columns\n across(\n # The cols argument takes a selection of columns to apply\n # a function to. Here, we want to apply the as.numeric()\n # function to all columns, so we use the function\n # everything() to select all columns.\n .cols = everything(),\n .fns = as.numeric\n )\n ) %>%\n # Convert the dataframe from wide to long format, so we have a\n # column for the time, a column for the state, and a column\n # for the proportion of the population in that state at that\n # time\n pivot_longer(\n # Don't pivot the time column\n cols = -c(time, iteration),\n names_to = \"state\",\n values_to = \"number\"\n ) %>%\n # Update the state column to be a factor, so the plot will\n # show the states in the correct order\n mutate(state = factor(state, levels = c(\"S\", \"I\", \"R\")))\n\nsir_colors <- c(S = \"#1f77b4\", I = \"#ff7f0e\", R = \"#FF3851\")\n\nggplot(\n gil_df,\n aes(\n x = time,\n y = number,\n color = state,\n group = interaction(iteration, state)\n )\n) +\n geom_line(linewidth = 1.5, alpha = .1) +\n scale_color_manual(values = sir_colors) +\n labs(\n x = \"Time\",\n y = \"Number\",\n color = \"State\"\n ) +\n guides(color = guide_legend(override.aes = list(alpha = 1))) +\n theme(legend.position = \"top\")\n\n\n\n\n\n\n\nThe Gillespie algorithm doesn’t cut any corners relative to the ODE model, but that comes at the cost of computational efficiency. For every event that happens (e.g. an infection) you also have to generate a random draw (e.g. a recovery) that you don’t use, except for comparison. And the bigger your population, the more possible events can happen. So the bigger the model (e.g. adding exposed or vaccinated classes, or heterogeneity, add transitions) and the bigger the population, the more calculation you need and and therefore the slower the model.", "crumbs": [ "Day 3", "11  R Session 04" diff --git a/_freeze/r-session-04/execute-results/html.json b/_freeze/r-session-04/execute-results/html.json index 7b48fbe..2cd3554 100644 --- a/_freeze/r-session-04/execute-results/html.json +++ b/_freeze/r-session-04/execute-results/html.json @@ -1,8 +1,8 @@ { - "hash": "f547f4be6e88b905aaecda04c0242ae8", + "hash": "b5991af5551cb01de54acc75877d60ce", "result": { "engine": "knitr", - "markdown": "---\nsubtitle: \"Stochastic SIR models\"\nabstract-title: \"\"\nabstract: |\n *review of basic algorithms for stochastic epidemic models*\nexecute:\n warning: false\nmetadata-files:\n - metadata/matthewferrari.yml\n - metadata/mathjax-packages.yml\neditor:\n markdown:\n wrap: sentence\n---\n\n# R Session 04\n\n:::{.callout-warning}\nThis session is under active development and is subject to change.\n:::\n\n## Setup\n\n\n::: {.cell}\n\n```{.r .cell-code}\nlibrary(here)\nlibrary(rio)\nlibrary(deSolve)\nlibrary(tidyverse)\n```\n:::\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\ntheme_set(theme_minimal())\n```\n:::\n\n\n:::callout-note\nThis R-session will go in parallel with the lecture on stochastic algorithms.\nFirst we'll go through the basics of each algorithm in the lecture and then we'll walk through the code and you can implement it for yourselves.\n:::\n\n## What is stochasticity and where does it come from?\n\nMuch of the world is uncertain (i.e. we don't know exactly how things work, what values are, or what tomorrow will hold).\nSome of that uncertainty is, at least theoretically, knowable and some is not.\nFor example, in our discussion of estimating $R_0$, there may be a very real $R_0$ for a given population and pathogen, even if we don't know it.\nThus, our estimate of $R_0$ may be \"uncertain\" (e.g. has a confidence interval around it, reflecting our certainty), but the models we've been developing so far are \"deterministic\", so conditional on a given value of $R_0$ the resulting epidemic curve is exactly specified by the model.\nIf we return to the code from R-session 1, we can plot a single deterministic realization of a model with $R_0 = 1.8$.\n\n\n::: {.cell}\n\n```{.r .cell-code}\nsir_model <- function(time, state, params, ...) {\n transmission <- params[\"transmission\"]\n recovery <- 1 / params[\"duration\"]\n\n S <- state[\"S\"]\n I <- state[\"I\"]\n R <- state[\"R\"]\n\n dSdt <- -transmission * S * I\n dIdt <- (transmission * S * I) - (recovery * I)\n dRdt <- recovery * I\n\n return(list(c(dSdt, dIdt, dRdt)))\n}\n```\n:::\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nsir_params <- c(transmission = 0.3, duration = 6)\nsir_init_states <- c(S = 0.99, I = 0.01, R = 0)\nsim_times <- seq(0, 200, by = 0.1)\n\nsir_sol <- ode(\n y = sir_init_states,\n times = sim_times,\n func = sir_model,\n parms = sir_params\n)\n```\n:::\n\n\n\n::: {.cell .column-body}\n\n```{.r .cell-code}\n# Turn the output from the ODE solver into a tibble (dataframe)\n# so we can manipulate and plot it easily\nsir_sol_df <- as_tibble(sir_sol) %>%\n # Convert all columns to numeric (they are currently type\n # deSolve so will produce warnings when plotting etc)\n mutate(\n # Rather than repeatedly type the same function for every\n # column, use the across() function to apply the function\n # to a selection of columns\n across(\n # The cols argument takes a selection of columns to apply\n # a function to. Here, we want to apply the as.numeric()\n # function to all columns, so we use the function\n # everything() to select all columns.\n .cols = everything(),\n .fns = as.numeric\n )\n ) %>%\n # Convert the dataframe from wide to long format, so we have a\n # column for the time, a column for the state, and a column\n # for the proportion of the population in that state at that\n # time\n pivot_longer(\n # Don't pivot the time column\n cols = -time,\n names_to = \"state\",\n values_to = \"proportion\"\n ) %>%\n # Update the state column to be a factor, so the plot will\n # show the states in the correct order\n mutate(state = factor(state, levels = c(\"S\", \"I\", \"R\")))\n\nsir_colors <- c(S = \"#1f77b4\", I = \"#ff7f0e\", R = \"#FF3851\")\n\nggplot(sir_sol_df, aes(x = time, y = proportion, color = state)) +\n geom_line(linewidth = 1.5) +\n scale_color_manual(values = sir_colors) +\n labs(\n x = \"Time\",\n y = \"Fraction\",\n color = \"State\"\n ) +\n theme(legend.position = \"top\")\n```\n\n::: {.cell-output-display}\n![](r-session-04_files/figure-html/unnamed-chunk-5-1.png){width=100%}\n:::\n:::\n\n\nNow, perhaps we used the tools from R-session 3 to estimate $R_0$ and we think a reasonable $95\\%$ confidence interval for $R_0$ is $(1.7,1.9)$.\nIf we're quite certain the duration of infection is 6 days, then that means our corresponding confidence interval on the transmission rate is $(.283,.317)$.\nAn entirely reasonable way to represent this uncertainty in the estimate of the transmission rate is to generate many random draws from within the confidence interval for the transmission rate and examine the resulting epidemic curves.\nFor this we can modify the code above with a loop that does this multiple times.\n\n\n::: {.cell}\n\n```{.r .cell-code}\n# the number of times we want to do this.\n# This is a matter of choice and computational capacity\n# (this model is small and quick, but that won't always be the case)\nnum_iterations <- 100\n\n# since we're going to do this num_iterations times we need a place\n# to store the results\n\nfor (iter in 1:num_iterations) {\n # each time take a random draw of the transmission rate from the\n # confidence interval\n sir_params <- c(transmission = runif(1, .283, .317), duration = 6)\n sir_init_states <- c(S = 0.99, I = 0.01, R = 0)\n sim_times <- seq(0, 200, by = 0.1)\n\n sir_sol <- ode(\n y = sir_init_states,\n times = sim_times,\n func = sir_model,\n parms = sir_params\n )\n sir_sol <- as_tibble(sir_sol)\n # mark this as the iter iteration of the loop\n sir_sol$iteration <- iter\n # if this is the first iteration, create a place to store the output\n # if it's the second or higher, append the output to the storage\n if (iter == 1) {\n sir_sol_storage <- sir_sol\n }\n if (iter > 1) {\n sir_sol_storage <- rbind(sir_sol_storage, sir_sol)\n }\n}\n```\n:::\n\n\n\n::: {.cell .column-body}\n\n```{.r .cell-code}\n# Turn the output from the ODE solver into a tibble (dataframe)\n# so we can manipulate and plot it easily\nsir_sol_df <- as_tibble(sir_sol_storage) %>%\n # Convert all columns to numeric (they are currently type deSolve\n # so will produce warnings when plotting etc)\n mutate(\n # Rather than repeatedly type the same function for every column,\n # use the across() function to apply the function to a selection\n # of columns\n across(\n # The cols argument takes a selection of columns to apply\n # a function to. Here, we want to apply the as.numeric()\n # function to all columns, so we use the function\n # everything() to select all columns.\n .cols = everything(),\n .fns = as.numeric\n )\n ) %>%\n # Convert the dataframe from wide to long format, so we have a\n # column for the time, a column for the state, and a column\n # for the proportion of the population in that state at that\n # time\n pivot_longer(\n # Don't pivot the time column\n cols = -c(time, iteration),\n names_to = \"state\",\n values_to = \"proportion\"\n ) %>%\n # Update the state column to be a factor, so the plot will\n # show the states in the correct order\n mutate(state = factor(state, levels = c(\"S\", \"I\", \"R\")))\n\nsir_colors <- c(S = \"#1f77b4\", I = \"#ff7f0e\", R = \"#FF3851\")\n\nggplot(\n sir_sol_df,\n aes(\n x = time,\n y = proportion,\n color = state,\n group = interaction(iteration, state)\n )\n) +\n geom_line(linewidth = 1.5, alpha = .1) +\n scale_color_manual(values = sir_colors) +\n labs(\n x = \"Time\",\n y = \"Fraction\",\n color = \"State\"\n ) +\n guides(color = guide_legend(override.aes = list(alpha = 1))) +\n theme(legend.position = \"top\")\n```\n\n::: {.cell-output-display}\n![](r-session-04_files/figure-html/unnamed-chunk-7-1.png){width=100%}\n:::\n:::\n\n\nHere note that we are still using our original code for the SIR model, so first we randomly draw the transmission rate, then conditional on that value, we run the deterministic model.\nNote that each run of the model still has the smooth curves, but each draw is a different set of smooth curves.\n\nThere are lots of ways you could extend this. Note that since we are drawing the transmission rate and fixing the duration of infection, the $R_0$ is slightly different in each simulation.\nYou could make each run have the same $R_0$ by recalculating the duration of infection based on the random draw of transmission rate; e.g. $R_0=1.8 = 0.29 * L$ means that $L=5.5$, so to ensure $R_0$ is the same in each run, you would need to change $L$ for each random draw.\nAlternatively, you might have uncertainty about \\emph{both} transmission and duration, so could for e.g. make random draws for both (which would again give a setting where $R_0$ varies from run to run).\nNone of these are more correct than another, the use case depends on which elements $R_0$, transmission, duration of infection, you are uncertain about.\n\nIn each of the above, the model is \\emph{deterministic}, meaning that for a given set of parameters, the resulting outbreak trajectory is always the same. For a \\emph{stochastic} model, each run of the model will vary, even if the parameters are the same.\nThis is because each event that occurs (e.g. someone getting infected or recovering) is analogous to a coin flip; even though the rules of the coin (the parameters) are always the same, the side it lands on is random.\nUnlike a deterministic model, to fully understand the behavior of a stochastic model, you need to run it more than one time, often many times, so these are necessarily more time consuming to work with.\n\nThere are many, many ways to make stochastic models and the steps can be way more complicated than flipping coins.\nBut there some foundational versions of these models that illustrate the trade-offs between exact representations of the random processes we think are happening and the computational time it takes to generate outputs. For what we'll do here, everything will be kind of fast, but in practice, for models of realistic scale, even a single stochastic run can take a while. And if you have to do thousands of runs, even shaving a few seconds or minutes can be important.\n\n## The Gillespie Algorithm\n\nThe Gillespie algorithm is the most explicit translation of the ODE form of the SIR model into a stochastic model.\nIt achieves this by noting that ODE-based models are written in terms of the rates at which events occur.\nAt any given point in time, the rate of all events in the ODE is known; what \\textbf{isn't} known, is which of the possible events will happen first.\nNote that even though we expect that the next thing to happen will be the thing that happens with the highest rate, it's possible that another thing may happen first by random chance.\nAnd then, since the rates in the SIR model are dependent on the value of the states (e.g. new infections depend on $\\beta$ and S and I) any change in the states then changes the rates and the likelihood of what will happen next.\n\nThe Gillespie algorithm proceeds by 1) calculating all the current rates, 2) randomly drawing exponential random variables that equate to the time until each event happens (recall we talked about the relationship between rate and time in the lectures), then 3) comparing those times and assuming that the next event to occur is the one that had the smallest randomly drawn time.\nThen you increment the states; e.g. an infection increases I by 1, reduces S by 1, and doesn't change R.\nImportantly, you then increment time forward by a step equal to the time until first event occurred.\nThen you recalculate the rates and randomly draw times, etc, and keep doing this over and over again.\nFor this SIR model, that keeps happening until you run out of infected individuals and there are no new events that can happen.\n\n\n::: {.cell}\n\n```{.r .cell-code}\n####################################################################\n# Parameters and initial conditions #\n####################################################################\n\nS <- 998 # number susceptible\nI <- 1 # number infected\nR <- 1 # number recovered\ntime <- 0\n\nbeta <- .5 # transmission rate\ngamma <- 1 / 7 # recovery rate\n```\n:::\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\n####################################################################\n# Gillespie Step #\n####################################################################\ngil_step <- function(SIR, beta, gamma) {\n # SIR is a vector containing 4 elements\n # S = scalar number of susceptibles now\n # I = scalar number of infecteds now\n # R = scalar number of recovereds now\n # time = current time\n # beta = transmission rate\n # gamma = recovery rate\n\n # draw two random exponential variables\n times <- rexp(\n 2,\n c(\n # the first is the rate of new infections: S*Beta*I/N\n beta * SIR[1] * SIR[2] / sum(SIR[1:3]),\n # the second is the rate of new recoveries: I*gamma\n SIR[2] * gamma\n )\n )\n\n return(list(\n # which.min identifies which of the two random variables is smallest,\n # and thus the first to happen: this is the transition that we'll make\n change_state = which.min(times),\n # this identifies how much time elapsed before the transition occurred,\n # so we can increment time forward\n time_step = min(times)\n ))\n}\n```\n:::\n\n\nThe gil_step() function increments the states forward by 1 event. Which event happens first, and the time it takes for that event to happen are both random variables. We then need to do this many times.\n\n\n::: {.cell}\n\n```{.r .cell-code}\n####################################################################\n# Simulate over time #\n####################################################################\n# here we set the random seed. This isn't necessary in general, but it allows us to write a \"random\" simulation\nset.seed(101)\ncounter <- 0 # set counter at 0\nwhile (all(I > 0)) {\n # continue until I is depleted\n counter <- counter + 1\n # counter for number of transitions: we don't know how many transitions\n # will happen until the simulation is over\n #\n # current SIR states\n sir_tmp <- c(S[counter], I[counter], R[counter], time[counter])\n step <- gil_step(sir_tmp, beta, gamma)\n if (step$change_state == 1) {\n # if transition is an infection, reduce S and increase I\n sir_tmp[1] <- sir_tmp[1] - 1\n sir_tmp[2] <- sir_tmp[2] + 1\n }\n if (step$change_state == 2) {\n # if transition is an recovery, reduce I and increase R\n sir_tmp[2] <- sir_tmp[2] - 1\n sir_tmp[3] <- sir_tmp[3] + 1\n }\n sir_tmp[4] <- sir_tmp[4] + step$time_step # increment time\n\n # Append changes\n S <- c(S, sir_tmp[1])\n I <- c(I, sir_tmp[2])\n R <- c(R, sir_tmp[3])\n time <- c(time, sir_tmp[4])\n\n # cat(S[counter],\\\"-\\\",I[counter],\\\"-\\\",R[counter],\\\"-\\\",time[counter], \\\".\\\\n\\\")\n # this prints the output as it goes\n # reset the seed so that every subsequent simulation IS random\n set.seed(NULL)\n}\n```\n:::\n\n\nWe can now plot the one realization of this stochastic outbreak. Notice that it has the same general shape as the deterministic outbreak, but is no longer smooth because each individual event over time happens randomly.\n\n\n::: {.cell .column-body}\n\n```{.r .cell-code}\n# Reuse the plotting code from above\n\n# Turn the output from gil_step() into a tibble (dataframe)\n# so we can manipulate and plot it easily\n# put elements in a data frame\ngil_df <- tibble(time, S, I, R) %>%\n # Convert all columns to numeric (they are currently type\n # deSolve so will produce warnings when plotting etc)\n mutate(\n # Rather than repeatedly type the same function for every\n # column, use the across() function to apply the function\n # to a selection of columns\n across(\n # The cols argument takes a selection of columns to apply\n # a function to. Here, we want to apply the as.numeric()\n # function to all columns, so we use the function\n # everything() to select all columns.\n .cols = everything(),\n .fns = as.numeric\n )\n ) %>%\n # Convert the dataframe from wide to long format, so we have a\n # column for the time, a column for the state, and a column\n # for the proportion of the population in that state at that\n # time\n pivot_longer(\n # Don't pivot the time column\n cols = -c(time),\n names_to = \"state\",\n values_to = \"number\"\n ) %>%\n # Update the state column to be a factor, so the plot will\n # show the states in the correct order\n mutate(state = factor(state, levels = c(\"S\", \"I\", \"R\")))\n\nsir_colors <- c(S = \"#1f77b4\", I = \"#ff7f0e\", R = \"#FF3851\")\n\nggplot(gil_df, aes(x = time, y = number, color = state)) +\n geom_line(linewidth = 1.5) +\n scale_color_manual(values = sir_colors) +\n labs(\n x = \"Time\",\n y = \"Number\",\n color = \"State\"\n ) +\n theme(legend.position = \"top\")\n```\n\n::: {.cell-output-display}\n![](r-session-04_files/figure-html/unnamed-chunk-11-1.png){width=100%}\n:::\n:::\n\n\nBecause each realization is stochastic, we need to generate many runs to see the general behavior.\nThe code below runs 10 iterations.\nGenerally, this would be considered a very small number of iterations.\nBut even for this very small model, running 100 or more means you'll be waiting for output.\nNote that this code is designed to be transparent not to be fast.\nMaking these run fast is beyond the scope of this assignment.\n\n\n::: {.cell}\n\n```{.r .cell-code}\n####################################################################\n# Simulate over time #\nnum_iterations <- 10 # number of realizations to simulate\n\nfor (iter in 1:num_iterations) {\n ####################################################################\n # Parameters and initial conditions #\n ####################################################################\n\n S <- 998 # number susceptible\n I <- 1 # number infected\n R <- 1 # number recovered\n time <- 0\n # initialize iteration counter\n iteration <- iter\n\n beta <- .5 # transmission rate\n gamma <- 1 / 7 # recovery rate\n\n ####################################################################\n # Simulate over time #\n ####################################################################\n counter <- 0 # set counter at 0\n while (all(I > 0)) {\n # continue until I is depleted\n counter <- counter + 1\n # counter for number of transitions: we don't know how many transitions\n # will happen until the simulation is over\n\n # current SIR states\n sir_tmp <- c(S[counter], I[counter], R[counter], time[counter])\n step <- gil_step(sir_tmp, beta, gamma)\n if (step$change_state == 1) {\n # if transition is an infection, reduce S and increase I\n sir_tmp[1] <- sir_tmp[1] - 1\n sir_tmp[2] <- sir_tmp[2] + 1\n }\n if (step$change_state == 2) {\n # if transition is an recovery, reduce I and increase R\n sir_tmp[2] <- sir_tmp[2] - 1\n sir_tmp[3] <- sir_tmp[3] + 1\n }\n # increment time\n sir_tmp[4] <- sir_tmp[4] + step$time_step\n\n # Append changes\n S <- c(S, sir_tmp[1])\n I <- c(I, sir_tmp[2])\n R <- c(R, sir_tmp[3])\n time <- c(time, sir_tmp[4])\n iteration <- c(iteration, iter)\n\n # cat(S[counter],\\\"-\\\",I[counter],\\\"-\\\",R[counter],\\\"-\\\",time[counter], \\\".\\\\n\\\")\n # this prints the output as it goes\n }\n if (iter == 1) {\n sir_gil_storage <- tibble(S, I, R, time, iteration)\n }\n if (iter > 1) {\n sir_gil_storage <- bind_rows(\n sir_gil_storage,\n tibble(S, I, R, time, iteration)\n )\n }\n}\n```\n:::\n\n\n\n::: {.cell .column-body}\n\n```{.r .cell-code}\n# Reuse the plotting code from above\n\n# Turn the output from gil_step() into a tibble (dataframe)\n# so we can manipulate and plot it easily\n# put elements in a data frame\ngil_df <- sir_gil_storage %>%\n # Convert all columns to numeric (they are currently type\n # deSolve so will produce warnings when plotting etc)\n mutate(\n # Rather than repeatedly type the same function for every\n # column, use the across() function to apply the function\n # to a selection of columns\n across(\n # The cols argument takes a selection of columns to apply\n # a function to. Here, we want to apply the as.numeric()\n # function to all columns, so we use the function\n # everything() to select all columns.\n .cols = everything(),\n .fns = as.numeric\n )\n ) %>%\n # Convert the dataframe from wide to long format, so we have a\n # column for the time, a column for the state, and a column\n # for the proportion of the population in that state at that\n # time\n pivot_longer(\n # Don't pivot the time column\n cols = -c(time, iteration),\n names_to = \"state\",\n values_to = \"number\"\n ) %>%\n # Update the state column to be a factor, so the plot will\n # show the states in the correct order\n mutate(state = factor(state, levels = c(\"S\", \"I\", \"R\")))\n\nsir_colors <- c(S = \"#1f77b4\", I = \"#ff7f0e\", R = \"#FF3851\")\n\nggplot(\n gil_df,\n aes(\n x = time,\n y = number,\n color = state,\n group = interaction(iteration, state)\n )\n) +\n geom_line(linewidth = 1.5, alpha = .1) +\n scale_color_manual(values = sir_colors) +\n labs(\n x = \"Time\",\n y = \"Number\",\n color = \"State\"\n ) +\n guides(color = guide_legend(override.aes = list(alpha = 1))) +\n theme(legend.position = \"top\")\n```\n\n::: {.cell-output-display}\n![](r-session-04_files/figure-html/unnamed-chunk-13-1.png){width=100%}\n:::\n:::\n\n\nThe Gillespie algorithm doesn't cut any corners relative to the ODE model, but that comes at the cost of computational efficiency.\nFor every event that happens (e.g. an infection) you also have to generate a random draw (e.g. a recovery) that you don't use, except for comparison.\nAnd the bigger your population, the more possible events can happen.\nSo the bigger the model (e.g. adding exposed or vaccinated classes, or heterogeneity, add transitions) and the bigger the population, the more calculation you need and and therefore the slower the model.\n\n\n## The Tau Leaping Algorithm\n\nThe Gillespie algorithm is great because it is an exact interpretation of the transitions in the ODE model: every change of state (e.g. infection or recovery) changes the rates for the next transition.\nDoing this comes at a computational cost.\nOne reasonable approximation is the Tau Leaping algorithm.\nHere, we move in discrete chunks of time and make random draws for multiple events occurring within that chunk.\nHere the computation scales with the number of time steps (and the number of states requiring transitions).\nBut, we rely on the result that, if rates stay constant, the number of events that occur in a discrete chunk of time can be approximated by a Poisson random variable.\nThus, we don't need to make random draws for every event.\nInstead we can make 1 draw for the multiple events that will occur in 1 day, or 1 week, etc.\nThe key here is the assumption that the rates stay constant; since each new infection or recovery will change the rates, we don't want \\emph{too many} to occur (in which case the rate at the start of the time step will be very different than the rate at the end of the time step).\nSo we are faced with a trade-off; very small time steps don't violate the assumptions, but the smaller the steps, the closer we are to Gillespie and the smaller the computational savings.\nThere's no right answer to what time step to use.\nHere we'll use 1 day, for convenience.\n\n\n::: {.cell}\n\n```{.r .cell-code}\n####################################################################\n# Parameters and initial conditions #\n####################################################################\n\nS <- 998 # number susceptible\nI <- 1 # number infected\nR <- 1 # number recovered\ntime <- 0\n\nbeta <- .5 # transmission rate\ngamma <- 1 / 7 # recovery rate\n```\n:::\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\n####################################################################\n# Tau Leaping Single time step #\n####################################################################\ntau_sir_step <- function(sims, S, I, R, beta, gamma, delta_t, ...) {\n # adapted from Aaron King's code\n # sims = number of simulations\n # S = initial susceptible population\n # I = initial infected population\n # R = initial recovered population\n # beta = transmission rate\n # gamma = recovery rate\n\n # total population size\n N <- S + I + R\n # new incident infections\n dSI <- rpois(n = sims, beta * S * (I / N) * delta_t)\n # recoveries\n dIR <- rpois(n = sims, gamma * I * delta_t)\n # note that this can be done with a binomial step as well\n # dSI <- rbinom(n=sims,size=S,prob=1-exp(-beta*(I/N)*delta_t))\n # new incident infections\n # dIR <- rbinom(n=sims,size=I,prob=1-exp(-gamma*delta_t))\n # recoveries\n\n # since it is possible for the transitions to drive the states negative,\n # we have to prevent that\n # change in S\n S <- pmax(S - dSI, 0)\n # change in I\n I <- pmax(I + dSI - dIR, 0)\n # change in R\n R <- R + dIR\n # note that dSI are the new incident infections\n cbind(S, I, R, dSI)\n}\n```\n:::\n\n\nNote in the code above that we're taking Poisson random draws, but there is some code commented out that uses binomial draws.\nThe former is exact for the theory, but it can give rise to settings where the transitions \"get ahead of themselves\" and you have more recoveries than infecteds, or more infections than susceptibles, which drives the states negative.\nWe can fix this by checking if the states go negative and disallowing this ... which is inelegant.\nWe can also fix this by using binomial draws, which are naturally constrained not to go negative.\nThe reason we don't automatically start with binomial draws is that they are computationally slower than Poisson draws (this occurs because the binomial distribution includes some combinatorial terms that are slow to compute).\nAs computers have gotten faster, this is less of an issue.\n\n\n::: {.cell}\n\n```{.r .cell-code}\n####################################################################\n# set up and storage for states #\n####################################################################\nmax_time <- 100\n# time to simulate over. With Tau Leaping the random draws scale with time\n# not with population size, so this is much more efficient than Gillespie\n# for large populations\nsims <- 1000\n# number of simulations: notice that we can do WAY more now\n\ns_mat <- matrix(S, 1, sims) # storage item for S for all simulations\ni_mat <- matrix(I, 1, sims) # storage item for I for all simulations\nr_mat <- matrix(R, 1, sims) # storage item for R for all simulations\nnew_cases <- matrix(0, 1, sims)\n# storage item for new cases (i.e. incidence) for all simulations\nn_mat <- S + I + R # storage item for N for all simulations\n```\n:::\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\n####################################################################\n# run over a time from 2 to T #\n####################################################################\n#\nfor (time_step in 2:max_time) {\n # loop over time, time_step is the index\n\n out <- tau_sir_step(\n sims,\n s_mat[time_step - 1, ],\n i_mat[time_step - 1, ],\n r_mat[time_step - 1, ],\n beta,\n gamma,\n delta_t = 1\n )\n # call to SIR step function above\n\n # update state\n s_mat <- rbind(s_mat, out[, 1])\n # update state\n i_mat <- rbind(i_mat, out[, 2])\n # update state\n r_mat <- rbind(r_mat, out[, 3])\n # update state -- note population size isn't changing, but this could be\n # updated with births/deaths\n n_mat <- rbind(n_mat, out[, 1] + out[, 2] + out[, 3])\n # update state\n new_cases <- rbind(new_cases, out[, 4])\n}\n```\n:::\n\n\nHopefully you can easily see that simulating 1000 iterations of Tau Leaping is \\emph{way} faster than Gillespie.\n\n\n::: {.cell .column-body}\n\n```{.r .cell-code}\n####################################################################\n# plotting #\n####################################################################\n# put output in a data frame\ntau_df <- tibble(\n S = array(s_mat),\n I = array(i_mat),\n R = array(r_mat),\n N = array(n_mat),\n cases = array(new_cases),\n time = rep(1:max_time, sims),\n iteration = rep(1:sims, each = max_time)\n)\n\ntau_df <- tau_df %>%\n pivot_longer(\n # Don't pivot the time column\n cols = -c(time, iteration),\n names_to = \"state\",\n values_to = \"number\"\n ) %>%\n mutate(state = factor(state, levels = c(\"S\", \"I\", \"R\", \"N\", \"cases\")))\n\n# fix this color\nsir_colors <- c(S = \"#1f77b4\", I = \"#ff7f0e\", R = \"#FF3851\", cases = \"#2ca02c\")\n\ntau_df %>%\n mutate(iteration = as.factor(iteration)) %>%\n filter(state %in% c(\"S\", \"I\", \"R\")) %>%\n ggplot(aes(\n x = time,\n y = number,\n group = interaction(iteration, state),\n color = state\n )) +\n geom_line(linewidth = 1.5, alpha = .1) +\n scale_color_manual(values = sir_colors) +\n labs(\n x = \"Time\",\n y = \"Number\",\n color = \"State\"\n ) +\n guides(color = guide_legend(override.aes = list(alpha = 1))) +\n theme(legend.position = \"top\")\n```\n\n::: {.cell-output-display}\n![](r-session-04_files/figure-html/unnamed-chunk-18-1.png){width=100%}\n:::\n:::\n\n\nAnd plotting the simulated trajectories should look pretty close to what we got with Gillespie.\nBut, because we can simulate many more iterations, we can start to observe some of the rarer behavior; e.g. even for simulations with $R_0 = 3.5$ there are some simulation runs for which the epidemic doesn't take off (S stays at 1000).\n\nNote that because we have stored the newly infected individuals as \"new cases\" then we can plot both the incidence (new cases) and prevalence (I) each day.\n\n\n::: {.cell .column-body}\n\n```{.r .cell-code}\n####################################################################\n# plotting #\n####################################################################\n# put output in a data frame\ntau_df <- tibble(\n S = array(s_mat),\n I = array(i_mat),\n R = array(r_mat),\n N = array(n_mat),\n cases = array(new_cases),\n time = rep(1:max_time, sims),\n iteration = rep(1:sims, each = max_time)\n)\n\ntau_df <- tau_df %>%\n pivot_longer(\n # Don't pivot the time column\n cols = -c(time, iteration),\n names_to = \"state\",\n values_to = \"number\"\n ) %>%\n mutate(state = factor(state, levels = c(\"S\", \"I\", \"R\", \"N\", \"cases\")))\n\n# fix this color\nsir_colors <- c(S = \"#1f77b4\", I = \"#ff7f0e\", R = \"#FF3851\", cases = \"#2ca02c\")\n\ntau_df %>%\n mutate(iteration = as.factor(iteration)) %>%\n filter(state %in% c(\"I\", \"cases\")) %>%\n ggplot(aes(\n x = time,\n y = number,\n group = interaction(iteration, state),\n color = state\n )) +\n geom_line(linewidth = 1.5, alpha = .1) +\n scale_color_manual(values = sir_colors) +\n labs(\n x = \"Time\",\n y = \"Number\",\n color = \"State\"\n ) +\n guides(color = guide_legend(override.aes = list(alpha = 1))) +\n theme(legend.position = \"top\")\n```\n\n::: {.cell-output-display}\n![](r-session-04_files/figure-html/unnamed-chunk-19-1.png){width=100%}\n:::\n:::\n\n\n:::callout-note\nRecall that the time series of new incident cases is very different than the time series of prevalent cases.\nRecall from our earlier discussion that the former are more likely to be what we would see in clinical surveillance.\nThe latter are what we might see if we did random testing in the population.\nWhich one would you expect to correspond best to environmental wastewater surveillance?\n::::\n\n## The Chain Binomial Algorithm\n\nIn the last section, we saw that by moving from continuous time to discrete time steps, we limit the complexity (the number of stochastic evaluations scales with time and the number of states, instead of population size and the number of states).\nYou may have also noticed that the individual steps in the Tau Leaping algorithm can be modeled as Poisson or Binomial random variables.\nCombining these two gives us another common algorithm for stochastic simulation, the Chain Binomial model.\nHere, we simplify further and use a time step that is equal to the infectious generation period; e.g. 2 weeks for measles, or 1 week for flu.\nIf we assume that time progresses in discrete, non-overlapping generations, then those that get infected at the start of one time step, recover at the end of that time step.\nThis simplification means that we only have to model the infection process as stochastic, and the recovery process at the end of the time step is now deterministic.\nThis is obviously unrealistic, but makes the model much simpler for formal statistical model fitting using likelihood and Bayesian methods.\nSo while it is imperfect, it remains as a common method for \"first-pass\" analyses.\n\nAs before, we start with initial conditions.\nHere we initialize with the same population as before, but notice that the transmission rate, $\\beta$ is different.\nTime is rescaled here to units of epidemic generation time, so $\\beta$ is scaled correspondingly and now, for this formulation, $\\beta$ is $R_0$.\n\n::: callout-note\nRecall that the basic reproduction number, $R_0$ is a function of the pathogen and the population of interest.\nIn models, both the population (e.g. mixing, heterogeneity, etc) and the pathogen (e.g. time-scale of recovery) are represented using approximations to the biology.\nHere, the choice of model forces a change in time-scale which means that we have to change $\\beta$ accordingly; thus, $\\beta$, which is ostensibly a biological parameter is also a property of the model.\n:::\n\n\n::: {.cell}\n\n```{.r .cell-code}\n####################################################################\n# Parameters and initial conditions #\n####################################################################\nS <- 998 # number susceptible\nI <- 1 # number infected\nR <- 1 # number recovered\ntime <- 0\n\nbeta <- 3.5 # transmission rate -- note the change in value\n```\n:::\n\n\nWhen we implement the forward simulation, we now only have to generate random draws for the infection process, which again simplifies the amount of stochastic simulation we have to do, and speeds up run times.\nWhile this might not be limiting in the scale of this activity, this DOES become an issue when we have to fit stochastic models.\nRecall from the lecture on parameter estimation that the basic recipe has us build a model and test out all, or many, values of the parameters (here $\\beta$) to find which one is closest to the data.\nHere, because each run of the stochastic simulation is different, we need to run simulations over many possible parameters, and for each parameter, we need to run many stochastic runs to characterize the average, or most likely, behavior.\nSo, the number of simulations can rapidly increase into the millions (bigger for models with more parameters), so every bit of savings in computational time can be valuable.\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\n####################################################################\n# Single time step #\n####################################################################\nchain_binomial_step <- function(sims, S, I, R, beta, ...) {\n # S = initial susceptible population\n # I = initial infected population\n # R = initial recovered population\n # beta is transmission rate\n\n # total population size\n N <- S + I + R\n # this is the only stochastic step; only do draw if I >0\n new_i <- rbinom(n = sims, pmax(S, 0), 1 - exp(-beta * pmax(I, 0) / N))\n # this step is now deterministic, because everyone from the past time\n # step recovers\n new_s <- S - new_i\n new_r <- R + I\n\n cbind(new_s, new_i, new_r)\n}\n```\n:::\n\n\nThen we can run the simulation over a set of discrete time steps, T.\nNote again that T is now the number of infectious generation times, since time is rescaled relative to the models above.\n\n\n::: {.cell}\n\n```{.r .cell-code}\n####################################################################\n# Run over many steps\t\t\t\t\t\t\t\t\t\t\t #\n####################################################################\n\nmax_time <- 20\n# note here that the time step is one infectious generation time,\n# so 7 days from gamma above\nsims <- 1000\n\ns_mat <- matrix(S, 1, sims) # storage item for S for all simulations\ni_mat <- matrix(I, 1, sims) # storage item for I for all simulations\nr_mat <- matrix(R, 1, sims) # storage item for R for all simulations\nn_mat <- S + I + R\n\nfor (time_step in 2:max_time) {\n out <- chain_binomial_step(\n sims,\n s_mat[time_step - 1, ],\n i_mat[time_step - 1, ],\n r_mat[time_step - 1, ],\n beta\n )\n # update state\n s_mat <- rbind(s_mat, out[, 1])\n # update state\n i_mat <- rbind(i_mat, out[, 2])\n # update state\n r_mat <- rbind(r_mat, out[, 3])\n # update state -- note population size isn't changing, but this could be\n # updated with births/deaths\n n_mat <- rbind(n_mat, out[, 1] + out[, 2] + out[, 3])\n}\n```\n:::\n\n\nWe can plot, but note again, that this is plotted in terms of epidemic generations on the X-axis, rather than days (as in the previous sections).\n\n\n::: {.cell .column-body}\n\n```{.r .cell-code}\n# put output in a data frame\nbin_df <- tibble(\n S = array(s_mat),\n I = array(i_mat),\n R = array(r_mat),\n N = array(matrix(n_mat)),\n time = rep(1:max_time, sims),\n iteration = rep(1:sims, each = max_time)\n)\n\nbin_df <- bin_df %>%\n pivot_longer(\n # Don't pivot the time column\n cols = -c(time, iteration),\n names_to = \"state\",\n values_to = \"number\"\n ) %>%\n mutate(state = factor(state, levels = c(\"S\", \"I\", \"R\", \"N\")))\n\n# fix this color\nsir_colors <- c(S = \"#1f77b4\", I = \"#ff7f0e\", R = \"#FF3851\")\n\nbin_df %>%\n mutate(iteration = as.factor(iteration)) %>%\n filter(state %in% c(\"S\", \"I\", \"R\")) %>%\n ggplot(aes(\n x = time,\n y = number,\n group = interaction(iteration, state),\n color = state\n )) +\n geom_line(linewidth = 1.5, alpha = .1) +\n scale_color_manual(values = sir_colors) +\n labs(\n x = \"Time\",\n y = \"Number\",\n color = \"State\"\n ) +\n guides(color = guide_legend(override.aes = list(alpha = 1))) +\n theme(legend.position = \"top\")\n```\n\n::: {.cell-output-display}\n![](r-session-04_files/figure-html/unnamed-chunk-23-1.png){width=100%}\n:::\n:::\n\n", + "markdown": "---\nsubtitle: \"Stochastic SIR models\"\nabstract-title: \"\"\nabstract: |\n *review of basic algorithms for stochastic epidemic models*\nexecute:\n warning: false\nmetadata-files:\n - metadata/matthewferrari.yml\n - metadata/mathjax-packages.yml\neditor:\n markdown:\n wrap: sentence\n---\n\n# R Session 04\n\n:::{.callout-warning}\nThis session is under active development and is subject to change.\n:::\n\n## Setup\n\n\n::: {.cell}\n\n```{.r .cell-code}\nlibrary(here)\nlibrary(rio)\nlibrary(deSolve)\nlibrary(tidyverse)\n```\n:::\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\ntheme_set(theme_minimal())\n```\n:::\n\n\n:::callout-note\nThis R-session will go in parallel with the lecture on stochastic algorithms.\nFirst we'll go through the basics of each algorithm in the lecture and then we'll walk through the code and you can implement it for yourselves.\n:::\n\n## What is stochasticity and where does it come from?\n\nMuch of the world is uncertain (i.e. we don't know exactly how things work, what values are, or what tomorrow will hold).\nSome of that uncertainty is, at least theoretically, knowable and some is not.\nFor example, in our discussion of estimating $R_0$, there may be a very real $R_0$ for a given population and pathogen, even if we don't know it.\nThus, our estimate of $R_0$ may be \"uncertain\" (e.g. has a confidence interval around it, reflecting our certainty), but the models we've been developing so far are \"deterministic\", so conditional on a given value of $R_0$ the resulting epidemic curve is exactly specified by the model.\nIf we return to the code from R-session 1, we can plot a single deterministic realization of a model with $R_0 = 1.8$.\n\n\n::: {.cell}\n\n```{.r .cell-code}\nsir_model <- function(time, state, params, ...) {\n transmission <- params[\"transmission\"]\n recovery <- 1 / params[\"duration\"]\n\n S <- state[\"S\"]\n I <- state[\"I\"]\n R <- state[\"R\"]\n\n dSdt <- -transmission * S * I\n dIdt <- (transmission * S * I) - (recovery * I)\n dRdt <- recovery * I\n\n return(list(c(dSdt, dIdt, dRdt)))\n}\n```\n:::\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nsir_params <- c(transmission = 0.3, duration = 6)\nsir_init_states <- c(S = 0.99, I = 0.01, R = 0)\nsim_times <- seq(0, 200, by = 0.1)\n\nsir_sol <- ode(\n y = sir_init_states,\n times = sim_times,\n func = sir_model,\n parms = sir_params\n)\n```\n:::\n\n\n\n::: {.cell .column-body}\n\n```{.r .cell-code}\n# Turn the output from the ODE solver into a tibble (dataframe)\n# so we can manipulate and plot it easily\nsir_sol_df <- as_tibble(sir_sol) %>%\n # Convert all columns to numeric (they are currently type\n # deSolve so will produce warnings when plotting etc)\n mutate(\n # Rather than repeatedly type the same function for every\n # column, use the across() function to apply the function\n # to a selection of columns\n across(\n # The cols argument takes a selection of columns to apply\n # a function to. Here, we want to apply the as.numeric()\n # function to all columns, so we use the function\n # everything() to select all columns.\n .cols = everything(),\n .fns = as.numeric\n )\n ) %>%\n # Convert the dataframe from wide to long format, so we have a\n # column for the time, a column for the state, and a column\n # for the proportion of the population in that state at that\n # time\n pivot_longer(\n # Don't pivot the time column\n cols = -time,\n names_to = \"state\",\n values_to = \"proportion\"\n ) %>%\n # Update the state column to be a factor, so the plot will\n # show the states in the correct order\n mutate(state = factor(state, levels = c(\"S\", \"I\", \"R\")))\n\nsir_colors <- c(S = \"#1f77b4\", I = \"#ff7f0e\", R = \"#FF3851\")\n\nggplot(sir_sol_df, aes(x = time, y = proportion, color = state)) +\n geom_line(linewidth = 1.5) +\n scale_color_manual(values = sir_colors) +\n labs(\n x = \"Time\",\n y = \"Fraction\",\n color = \"State\"\n ) +\n theme(legend.position = \"top\")\n```\n\n::: {.cell-output-display}\n![](r-session-04_files/figure-html/unnamed-chunk-5-1.png){width=100%}\n:::\n:::\n\n\nNow, perhaps we used the tools from R-session 3 to estimate $R_0$ and we think a reasonable $95\\%$ confidence interval for $R_0$ is $(1.7,1.9)$.\nIf we're quite certain the duration of infection is 6 days, then that means our corresponding confidence interval on the transmission rate is $(.283,.317)$.\nAn entirely reasonable way to represent this uncertainty in the estimate of the transmission rate is to generate many random draws from within the confidence interval for the transmission rate and examine the resulting epidemic curves.\nFor this we can modify the code above with a loop that does this multiple times.\n\n\n::: {.cell}\n\n```{.r .cell-code}\n# the number of times we want to do this.\n# This is a matter of choice and computational capacity\n# (this model is small and quick, but that won't always be the case)\nnum_iterations <- 100\n\n# since we're going to do this num_iterations times we need a place\n# to store the results\n\nfor (iter in 1:num_iterations) {\n # each time take a random draw of the transmission rate from the\n # confidence interval\n sir_params <- c(transmission = runif(1, .283, .317), duration = 6)\n sir_init_states <- c(S = 0.99, I = 0.01, R = 0)\n sim_times <- seq(0, 200, by = 0.1)\n\n sir_sol <- ode(\n y = sir_init_states,\n times = sim_times,\n func = sir_model,\n parms = sir_params\n )\n sir_sol <- as_tibble(sir_sol)\n # mark this as the iter iteration of the loop\n sir_sol$iteration <- iter\n # if this is the first iteration, create a place to store the output\n # if it's the second or higher, append the output to the storage\n if (iter == 1) {\n sir_sol_storage <- sir_sol\n }\n if (iter > 1) {\n sir_sol_storage <- rbind(sir_sol_storage, sir_sol)\n }\n}\n```\n:::\n\n\n\n::: {.cell .column-body}\n\n```{.r .cell-code}\n# Turn the output from the ODE solver into a tibble (dataframe)\n# so we can manipulate and plot it easily\nsir_sol_df <- as_tibble(sir_sol_storage) %>%\n # Convert all columns to numeric (they are currently type deSolve\n # so will produce warnings when plotting etc)\n mutate(\n # Rather than repeatedly type the same function for every column,\n # use the across() function to apply the function to a selection\n # of columns\n across(\n # The cols argument takes a selection of columns to apply\n # a function to. Here, we want to apply the as.numeric()\n # function to all columns, so we use the function\n # everything() to select all columns.\n .cols = everything(),\n .fns = as.numeric\n )\n ) %>%\n # Convert the dataframe from wide to long format, so we have a\n # column for the time, a column for the state, and a column\n # for the proportion of the population in that state at that\n # time\n pivot_longer(\n # Don't pivot the time column\n cols = -c(time, iteration),\n names_to = \"state\",\n values_to = \"proportion\"\n ) %>%\n # Update the state column to be a factor, so the plot will\n # show the states in the correct order\n mutate(state = factor(state, levels = c(\"S\", \"I\", \"R\")))\n\nsir_colors <- c(S = \"#1f77b4\", I = \"#ff7f0e\", R = \"#FF3851\")\n\nggplot(\n sir_sol_df,\n aes(\n x = time,\n y = proportion,\n color = state,\n group = interaction(iteration, state)\n )\n) +\n geom_line(linewidth = 1.5, alpha = .1) +\n scale_color_manual(values = sir_colors) +\n labs(\n x = \"Time\",\n y = \"Fraction\",\n color = \"State\"\n ) +\n guides(color = guide_legend(override.aes = list(alpha = 1))) +\n theme(legend.position = \"top\")\n```\n\n::: {.cell-output-display}\n![](r-session-04_files/figure-html/unnamed-chunk-7-1.png){width=100%}\n:::\n:::\n\n\nHere note that we are still using our original code for the SIR model, so first we randomly draw the transmission rate, then conditional on that value, we run the deterministic model.\nNote that each run of the model still has the smooth curves, but each draw is a different set of smooth curves.\n\nThere are lots of ways you could extend this. Note that since we are drawing the transmission rate and fixing the duration of infection, the $R_0$ is slightly different in each simulation.\nYou could make each run have the same $R_0$ by recalculating the duration of infection based on the random draw of transmission rate; e.g. $R_0=1.8 = 0.29 * L$ means that $L=5.5$, so to ensure $R_0$ is the same in each run, you would need to change $L$ for each random draw.\nAlternatively, you might have uncertainty about \\emph{both} transmission and duration, so could for e.g. make random draws for both (which would again give a setting where $R_0$ varies from run to run).\nNone of these are more correct than another, the use case depends on which elements $R_0$, transmission, duration of infection, you are uncertain about.\n\nIn each of the above, the model is \\emph{deterministic}, meaning that for a given set of parameters, the resulting outbreak trajectory is always the same. For a \\emph{stochastic} model, each run of the model will vary, even if the parameters are the same.\nThis is because each event that occurs (e.g. someone getting infected or recovering) is analogous to a coin flip; even though the rules of the coin (the parameters) are always the same, the side it lands on is random.\nUnlike a deterministic model, to fully understand the behavior of a stochastic model, you need to run it more than one time, often many times, so these are necessarily more time consuming to work with.\n\nThere are many, many ways to make stochastic models and the steps can be way more complicated than flipping coins.\nBut there some foundational versions of these models that illustrate the trade-offs between exact representations of the random processes we think are happening and the computational time it takes to generate outputs. For what we'll do here, everything will be kind of fast, but in practice, for models of realistic scale, even a single stochastic run can take a while. And if you have to do thousands of runs, even shaving a few seconds or minutes can be important.\n\n## The Gillespie Algorithm\n\nThe Gillespie algorithm is the most explicit translation of the ODE form of the SIR model into a stochastic model.\nIt achieves this by noting that ODE-based models are written in terms of the rates at which events occur.\nAt any given point in time, the rate of all events in the ODE is known; what \\textbf{isn't} known, is which of the possible events will happen first.\nNote that even though we expect that the next thing to happen will be the thing that happens with the highest rate, it's possible that another thing may happen first by random chance.\nAnd then, since the rates in the SIR model are dependent on the value of the states (e.g. new infections depend on $\\beta$ and S and I) any change in the states then changes the rates and the likelihood of what will happen next.\n\nThe Gillespie algorithm proceeds by 1) calculating all the current rates, 2) randomly drawing exponential random variables that equate to the time until each event happens (recall we talked about the relationship between rate and time in the lectures), then 3) comparing those times and assuming that the next event to occur is the one that had the smallest randomly drawn time.\nThen you increment the states; e.g. an infection increases I by 1, reduces S by 1, and doesn't change R.\nImportantly, you then increment time forward by a step equal to the time until first event occurred.\nThen you recalculate the rates and randomly draw times, etc, and keep doing this over and over again.\nFor this SIR model, that keeps happening until you run out of infected individuals and there are no new events that can happen.\n\n\n::: {.cell}\n\n```{.r .cell-code}\n####################################################################\n# Parameters and initial conditions #\n####################################################################\n\nS <- 998 # number susceptible\nI <- 1 # number infected\nR <- 1 # number recovered\ntime <- 0\n\nbeta <- 0.5 # transmission rate\ngamma <- 1 / 7 # recovery rate\n```\n:::\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\n####################################################################\n# Gillespie Step #\n####################################################################\ngil_step <- function(SIR, beta, gamma) {\n # SIR is a vector containing 4 elements\n # S = scalar number of susceptibles now\n # I = scalar number of infecteds now\n # R = scalar number of recovereds now\n # time = current time\n # beta = transmission rate\n # gamma = recovery rate\n\n # draw two random exponential variables\n times <- rexp(\n 2,\n c(\n # the first is the rate of new infections: S*Beta*I/N\n beta * SIR[1] * SIR[2] / sum(SIR[1:3]),\n # the second is the rate of new recoveries: I*gamma\n SIR[2] * gamma\n )\n )\n\n return(list(\n # which.min identifies which of the two random variables is smallest,\n # and thus the first to happen: this is the transition that we'll make\n change_state = which.min(times),\n # this identifies how much time elapsed before the transition occurred,\n # so we can increment time forward\n time_step = min(times)\n ))\n}\n```\n:::\n\n\nThe gil_step() function increments the states forward by 1 event. Which event happens first, and the time it takes for that event to happen are both random variables. We then need to do this many times.\n\n\n::: {.cell}\n\n```{.r .cell-code}\n####################################################################\n# Simulate over time #\n####################################################################\n# here we set the random seed. This isn't necessary in general, but it allows us to write a \"random\" simulation\nset.seed(101)\ncounter <- 0 # set counter at 0\nwhile (all(I > 0)) {\n # continue until I is depleted\n counter <- counter + 1\n # counter for number of transitions: we don't know how many transitions\n # will happen until the simulation is over\n #\n # current SIR states\n sir_tmp <- c(S[counter], I[counter], R[counter], time[counter])\n step <- gil_step(sir_tmp, beta, gamma)\n if (step$change_state == 1) {\n # if transition is an infection, reduce S and increase I\n sir_tmp[1] <- sir_tmp[1] - 1\n sir_tmp[2] <- sir_tmp[2] + 1\n }\n if (step$change_state == 2) {\n # if transition is an recovery, reduce I and increase R\n sir_tmp[2] <- sir_tmp[2] - 1\n sir_tmp[3] <- sir_tmp[3] + 1\n }\n sir_tmp[4] <- sir_tmp[4] + step$time_step # increment time\n\n # Append changes\n S <- c(S, sir_tmp[1])\n I <- c(I, sir_tmp[2])\n R <- c(R, sir_tmp[3])\n time <- c(time, sir_tmp[4])\n\n # cat(S[counter],\\\"-\\\",I[counter],\\\"-\\\",R[counter],\\\"-\\\",time[counter], \\\".\\\\n\\\")\n # this prints the output as it goes\n # reset the seed so that every subsequent simulation IS random\n set.seed(NULL)\n}\n```\n:::\n\n\nWe can now plot the one realization of this stochastic outbreak. Notice that it has the same general shape as the deterministic outbreak, but is no longer smooth because each individual event over time happens randomly.\n\n\n::: {.cell .column-body}\n\n```{.r .cell-code}\n# Reuse the plotting code from above\n\n# Turn the output from gil_step() into a tibble (dataframe)\n# so we can manipulate and plot it easily\n# put elements in a data frame\ngil_df <- tibble(time, S, I, R) %>%\n # Convert all columns to numeric (they are currently type\n # deSolve so will produce warnings when plotting etc)\n mutate(\n # Rather than repeatedly type the same function for every\n # column, use the across() function to apply the function\n # to a selection of columns\n across(\n # The cols argument takes a selection of columns to apply\n # a function to. Here, we want to apply the as.numeric()\n # function to all columns, so we use the function\n # everything() to select all columns.\n .cols = everything(),\n .fns = as.numeric\n )\n ) %>%\n # Convert the dataframe from wide to long format, so we have a\n # column for the time, a column for the state, and a column\n # for the proportion of the population in that state at that\n # time\n pivot_longer(\n # Don't pivot the time column\n cols = -c(time),\n names_to = \"state\",\n values_to = \"number\"\n ) %>%\n # Update the state column to be a factor, so the plot will\n # show the states in the correct order\n mutate(state = factor(state, levels = c(\"S\", \"I\", \"R\")))\n\nsir_colors <- c(S = \"#1f77b4\", I = \"#ff7f0e\", R = \"#FF3851\")\n\nggplot(gil_df, aes(x = time, y = number, color = state)) +\n geom_line(linewidth = 1.5) +\n scale_color_manual(values = sir_colors) +\n labs(\n x = \"Time\",\n y = \"Number\",\n color = \"State\"\n ) +\n theme(legend.position = \"top\")\n```\n\n::: {.cell-output-display}\n![](r-session-04_files/figure-html/unnamed-chunk-11-1.png){width=100%}\n:::\n:::\n\n\nBecause each realization is stochastic, we need to generate many runs to see the general behavior.\nThe code below runs 10 iterations.\nGenerally, this would be considered a very small number of iterations.\nBut even for this very small model, running 100 or more means you'll be waiting for output.\nNote that this code is designed to be transparent not to be fast.\nMaking these run fast is beyond the scope of this assignment.\n\n\n::: {.cell}\n\n```{.r .cell-code}\n####################################################################\n# Simulate over time #\nnum_iterations <- 10 # number of realizations to simulate\n\nfor (iter in 1:num_iterations) {\n ####################################################################\n # Parameters and initial conditions #\n ####################################################################\n\n S <- 998 # number susceptible\n I <- 1 # number infected\n R <- 1 # number recovered\n time <- 0\n # initialize iteration counter\n iteration <- iter\n\n beta <- .5 # transmission rate\n gamma <- 1 / 7 # recovery rate\n\n ####################################################################\n # Simulate over time #\n ####################################################################\n counter <- 0 # set counter at 0\n while (all(I > 0)) {\n # continue until I is depleted\n counter <- counter + 1\n # counter for number of transitions: we don't know how many transitions\n # will happen until the simulation is over\n\n # current SIR states\n sir_tmp <- c(S[counter], I[counter], R[counter], time[counter])\n step <- gil_step(sir_tmp, beta, gamma)\n if (step$change_state == 1) {\n # if transition is an infection, reduce S and increase I\n sir_tmp[1] <- sir_tmp[1] - 1\n sir_tmp[2] <- sir_tmp[2] + 1\n }\n if (step$change_state == 2) {\n # if transition is an recovery, reduce I and increase R\n sir_tmp[2] <- sir_tmp[2] - 1\n sir_tmp[3] <- sir_tmp[3] + 1\n }\n # increment time\n sir_tmp[4] <- sir_tmp[4] + step$time_step\n\n # Append changes\n S <- c(S, sir_tmp[1])\n I <- c(I, sir_tmp[2])\n R <- c(R, sir_tmp[3])\n time <- c(time, sir_tmp[4])\n iteration <- c(iteration, iter)\n\n # cat(S[counter],\\\"-\\\",I[counter],\\\"-\\\",R[counter],\\\"-\\\",time[counter], \\\".\\\\n\\\")\n # this prints the output as it goes\n }\n if (iter == 1) {\n sir_gil_storage <- tibble(S, I, R, time, iteration)\n }\n if (iter > 1) {\n sir_gil_storage <- bind_rows(\n sir_gil_storage,\n tibble(S, I, R, time, iteration)\n )\n }\n}\n```\n:::\n\n\n\n::: {.cell .column-body}\n\n```{.r .cell-code}\n# Reuse the plotting code from above\n\n# Turn the output from gil_step() into a tibble (dataframe)\n# so we can manipulate and plot it easily\n# put elements in a data frame\ngil_df <- sir_gil_storage %>%\n # Convert all columns to numeric (they are currently type\n # deSolve so will produce warnings when plotting etc)\n mutate(\n # Rather than repeatedly type the same function for every\n # column, use the across() function to apply the function\n # to a selection of columns\n across(\n # The cols argument takes a selection of columns to apply\n # a function to. Here, we want to apply the as.numeric()\n # function to all columns, so we use the function\n # everything() to select all columns.\n .cols = everything(),\n .fns = as.numeric\n )\n ) %>%\n # Convert the dataframe from wide to long format, so we have a\n # column for the time, a column for the state, and a column\n # for the proportion of the population in that state at that\n # time\n pivot_longer(\n # Don't pivot the time column\n cols = -c(time, iteration),\n names_to = \"state\",\n values_to = \"number\"\n ) %>%\n # Update the state column to be a factor, so the plot will\n # show the states in the correct order\n mutate(state = factor(state, levels = c(\"S\", \"I\", \"R\")))\n\nsir_colors <- c(S = \"#1f77b4\", I = \"#ff7f0e\", R = \"#FF3851\")\n\nggplot(\n gil_df,\n aes(\n x = time,\n y = number,\n color = state,\n group = interaction(iteration, state)\n )\n) +\n geom_line(linewidth = 1.5, alpha = .1) +\n scale_color_manual(values = sir_colors) +\n labs(\n x = \"Time\",\n y = \"Number\",\n color = \"State\"\n ) +\n guides(color = guide_legend(override.aes = list(alpha = 1))) +\n theme(legend.position = \"top\")\n```\n\n::: {.cell-output-display}\n![](r-session-04_files/figure-html/unnamed-chunk-13-1.png){width=100%}\n:::\n:::\n\n\nThe Gillespie algorithm doesn't cut any corners relative to the ODE model, but that comes at the cost of computational efficiency.\nFor every event that happens (e.g. an infection) you also have to generate a random draw (e.g. a recovery) that you don't use, except for comparison.\nAnd the bigger your population, the more possible events can happen.\nSo the bigger the model (e.g. adding exposed or vaccinated classes, or heterogeneity, add transitions) and the bigger the population, the more calculation you need and and therefore the slower the model.\n\n\n## The Tau Leaping Algorithm\n\nThe Gillespie algorithm is great because it is an exact interpretation of the transitions in the ODE model: every change of state (e.g. infection or recovery) changes the rates for the next transition.\nDoing this comes at a computational cost.\nOne reasonable approximation is the Tau Leaping algorithm.\nHere, we move in discrete chunks of time and make random draws for multiple events occurring within that chunk.\nHere the computation scales with the number of time steps (and the number of states requiring transitions).\nBut, we rely on the result that, if rates stay constant, the number of events that occur in a discrete chunk of time can be approximated by a Poisson random variable.\nThus, we don't need to make random draws for every event.\nInstead we can make 1 draw for the multiple events that will occur in 1 day, or 1 week, etc.\nThe key here is the assumption that the rates stay constant; since each new infection or recovery will change the rates, we don't want \\emph{too many} to occur (in which case the rate at the start of the time step will be very different than the rate at the end of the time step).\nSo we are faced with a trade-off; very small time steps don't violate the assumptions, but the smaller the steps, the closer we are to Gillespie and the smaller the computational savings.\nThere's no right answer to what time step to use.\nHere we'll use 1 day, for convenience.\n\n\n::: {.cell}\n\n```{.r .cell-code}\n####################################################################\n# Parameters and initial conditions #\n####################################################################\n\nS <- 998 # number susceptible\nI <- 1 # number infected\nR <- 1 # number recovered\ntime <- 0\n\nbeta <- .5 # transmission rate\ngamma <- 1 / 7 # recovery rate\n```\n:::\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\n####################################################################\n# Tau Leaping Single time step #\n####################################################################\ntau_sir_step <- function(sims, S, I, R, beta, gamma, delta_t, ...) {\n # adapted from Aaron King's code\n # sims = number of simulations\n # S = initial susceptible population\n # I = initial infected population\n # R = initial recovered population\n # beta = transmission rate\n # gamma = recovery rate\n\n # total population size\n N <- S + I + R\n # new incident infections\n dSI <- rpois(n = sims, beta * S * (I / N) * delta_t)\n # recoveries\n dIR <- rpois(n = sims, gamma * I * delta_t)\n # note that this can be done with a binomial step as well\n # dSI <- rbinom(n=sims,size=S,prob=1-exp(-beta*(I/N)*delta_t))\n # new incident infections\n # dIR <- rbinom(n=sims,size=I,prob=1-exp(-gamma*delta_t))\n # recoveries\n\n # since it is possible for the transitions to drive the states negative,\n # we have to prevent that\n # change in S\n S <- pmax(S - dSI, 0)\n # change in I\n I <- pmax(I + dSI - dIR, 0)\n # change in R\n R <- R + dIR\n # note that dSI are the new incident infections\n cbind(S, I, R, dSI)\n}\n```\n:::\n\n\nNote in the code above that we're taking Poisson random draws, but there is some code commented out that uses binomial draws.\nThe former is exact for the theory, but it can give rise to settings where the transitions \"get ahead of themselves\" and you have more recoveries than infecteds, or more infections than susceptibles, which drives the states negative.\nWe can fix this by checking if the states go negative and disallowing this ... which is inelegant.\nWe can also fix this by using binomial draws, which are naturally constrained not to go negative.\nThe reason we don't automatically start with binomial draws is that they are computationally slower than Poisson draws (this occurs because the binomial distribution includes some combinatorial terms that are slow to compute).\nAs computers have gotten faster, this is less of an issue.\n\n\n::: {.cell}\n\n```{.r .cell-code}\n####################################################################\n# set up and storage for states #\n####################################################################\nmax_time <- 100\n# time to simulate over. With Tau Leaping the random draws scale with time\n# not with population size, so this is much more efficient than Gillespie\n# for large populations\nsims <- 1000\n# number of simulations: notice that we can do WAY more now\n\ns_mat <- matrix(S, 1, sims) # storage item for S for all simulations\ni_mat <- matrix(I, 1, sims) # storage item for I for all simulations\nr_mat <- matrix(R, 1, sims) # storage item for R for all simulations\nnew_cases <- matrix(0, 1, sims)\n# storage item for new cases (i.e. incidence) for all simulations\nn_mat <- S + I + R # storage item for N for all simulations\n```\n:::\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\n####################################################################\n# run over a time from 2 to T #\n####################################################################\n#\nfor (time_step in 2:max_time) {\n # loop over time, time_step is the index\n\n out <- tau_sir_step(\n sims,\n s_mat[time_step - 1, ],\n i_mat[time_step - 1, ],\n r_mat[time_step - 1, ],\n beta,\n gamma,\n delta_t = 1\n )\n # call to SIR step function above\n\n # update state\n s_mat <- rbind(s_mat, out[, 1])\n # update state\n i_mat <- rbind(i_mat, out[, 2])\n # update state\n r_mat <- rbind(r_mat, out[, 3])\n # update state -- note population size isn't changing, but this could be\n # updated with births/deaths\n n_mat <- rbind(n_mat, out[, 1] + out[, 2] + out[, 3])\n # update state\n new_cases <- rbind(new_cases, out[, 4])\n}\n```\n:::\n\n\nHopefully you can easily see that simulating 1000 iterations of Tau Leaping is \\emph{way} faster than Gillespie.\n\n\n::: {.cell .column-body}\n\n```{.r .cell-code}\n####################################################################\n# plotting #\n####################################################################\n# put output in a data frame\ntau_df <- tibble(\n S = array(s_mat),\n I = array(i_mat),\n R = array(r_mat),\n N = array(n_mat),\n cases = array(new_cases),\n time = rep(1:max_time, sims),\n iteration = rep(1:sims, each = max_time)\n)\n\ntau_df <- tau_df %>%\n pivot_longer(\n # Don't pivot the time column\n cols = -c(time, iteration),\n names_to = \"state\",\n values_to = \"number\"\n ) %>%\n mutate(state = factor(state, levels = c(\"S\", \"I\", \"R\", \"N\", \"cases\")))\n\n# fix this color\nsir_colors <- c(S = \"#1f77b4\", I = \"#ff7f0e\", R = \"#FF3851\", cases = \"#2ca02c\")\n\ntau_df %>%\n mutate(iteration = as.factor(iteration)) %>%\n filter(state %in% c(\"S\", \"I\", \"R\")) %>%\n ggplot(aes(\n x = time,\n y = number,\n group = interaction(iteration, state),\n color = state\n )) +\n geom_line(linewidth = 1.5, alpha = .1) +\n scale_color_manual(values = sir_colors) +\n labs(\n x = \"Time\",\n y = \"Number\",\n color = \"State\"\n ) +\n guides(color = guide_legend(override.aes = list(alpha = 1))) +\n theme(legend.position = \"top\")\n```\n\n::: {.cell-output-display}\n![](r-session-04_files/figure-html/unnamed-chunk-18-1.png){width=100%}\n:::\n:::\n\n\nAnd plotting the simulated trajectories should look pretty close to what we got with Gillespie.\nBut, because we can simulate many more iterations, we can start to observe some of the rarer behavior; e.g. even for simulations with $R_0 = 3.5$ there are some simulation runs for which the epidemic doesn't take off (S stays at 1000).\n\nNote that because we have stored the newly infected individuals as \"new cases\" then we can plot both the incidence (new cases) and prevalence (I) each day.\n\n\n::: {.cell .column-body}\n\n```{.r .cell-code}\n####################################################################\n# plotting #\n####################################################################\n# put output in a data frame\ntau_df <- tibble(\n S = array(s_mat),\n I = array(i_mat),\n R = array(r_mat),\n N = array(n_mat),\n cases = array(new_cases),\n time = rep(1:max_time, sims),\n iteration = rep(1:sims, each = max_time)\n)\n\ntau_df <- tau_df %>%\n pivot_longer(\n # Don't pivot the time column\n cols = -c(time, iteration),\n names_to = \"state\",\n values_to = \"number\"\n ) %>%\n mutate(state = factor(state, levels = c(\"S\", \"I\", \"R\", \"N\", \"cases\")))\n\n# fix this color\nsir_colors <- c(S = \"#1f77b4\", I = \"#ff7f0e\", R = \"#FF3851\", cases = \"#2ca02c\")\n\ntau_df %>%\n mutate(iteration = as.factor(iteration)) %>%\n filter(state %in% c(\"I\", \"cases\")) %>%\n ggplot(aes(\n x = time,\n y = number,\n group = interaction(iteration, state),\n color = state\n )) +\n geom_line(linewidth = 1.5, alpha = .1) +\n scale_color_manual(values = sir_colors) +\n labs(\n x = \"Time\",\n y = \"Number\",\n color = \"State\"\n ) +\n guides(color = guide_legend(override.aes = list(alpha = 1))) +\n theme(legend.position = \"top\")\n```\n\n::: {.cell-output-display}\n![](r-session-04_files/figure-html/unnamed-chunk-19-1.png){width=100%}\n:::\n:::\n\n\n:::callout-note\nRecall that the time series of new incident cases is very different than the time series of prevalent cases.\nRecall from our earlier discussion that the former are more likely to be what we would see in clinical surveillance.\nThe latter are what we might see if we did random testing in the population.\nWhich one would you expect to correspond best to environmental wastewater surveillance?\n::::\n\n## The Chain Binomial Algorithm\n\nIn the last section, we saw that by moving from continuous time to discrete time steps, we limit the complexity (the number of stochastic evaluations scales with time and the number of states, instead of population size and the number of states).\nYou may have also noticed that the individual steps in the Tau Leaping algorithm can be modeled as Poisson or Binomial random variables.\nCombining these two gives us another common algorithm for stochastic simulation, the Chain Binomial model.\nHere, we simplify further and use a time step that is equal to the infectious generation period; e.g. 2 weeks for measles, or 1 week for flu.\nIf we assume that time progresses in discrete, non-overlapping generations, then those that get infected at the start of one time step, recover at the end of that time step.\nThis simplification means that we only have to model the infection process as stochastic, and the recovery process at the end of the time step is now deterministic.\nThis is obviously unrealistic, but makes the model much simpler for formal statistical model fitting using likelihood and Bayesian methods.\nSo while it is imperfect, it remains as a common method for \"first-pass\" analyses.\n\nAs before, we start with initial conditions.\nHere we initialize with the same population as before, but notice that the transmission rate, $\\beta$ is different.\nTime is rescaled here to units of epidemic generation time, so $\\beta$ is scaled correspondingly and now, for this formulation, $\\beta$ is $R_0$.\n\n::: callout-note\nRecall that the basic reproduction number, $R_0$ is a function of the pathogen and the population of interest.\nIn models, both the population (e.g. mixing, heterogeneity, etc) and the pathogen (e.g. time-scale of recovery) are represented using approximations to the biology.\nHere, the choice of model forces a change in time-scale which means that we have to change $\\beta$ accordingly; thus, $\\beta$, which is ostensibly a biological parameter is also a property of the model.\n:::\n\n\n::: {.cell}\n\n```{.r .cell-code}\n####################################################################\n# Parameters and initial conditions #\n####################################################################\nS <- 998 # number susceptible\nI <- 1 # number infected\nR <- 1 # number recovered\ntime <- 0\n\nbeta <- 3.5 # transmission rate -- note the change in value\n```\n:::\n\n\nWhen we implement the forward simulation, we now only have to generate random draws for the infection process, which again simplifies the amount of stochastic simulation we have to do, and speeds up run times.\nWhile this might not be limiting in the scale of this activity, this DOES become an issue when we have to fit stochastic models.\nRecall from the lecture on parameter estimation that the basic recipe has us build a model and test out all, or many, values of the parameters (here $\\beta$) to find which one is closest to the data.\nHere, because each run of the stochastic simulation is different, we need to run simulations over many possible parameters, and for each parameter, we need to run many stochastic runs to characterize the average, or most likely, behavior.\nSo, the number of simulations can rapidly increase into the millions (bigger for models with more parameters), so every bit of savings in computational time can be valuable.\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\n####################################################################\n# Single time step #\n####################################################################\nchain_binomial_step <- function(sims, S, I, R, beta, ...) {\n # S = initial susceptible population\n # I = initial infected population\n # R = initial recovered population\n # beta is transmission rate\n\n # total population size\n N <- S + I + R\n # this is the only stochastic step; only do draw if I >0\n new_i <- rbinom(n = sims, pmax(S, 0), 1 - exp(-beta * pmax(I, 0) / N))\n # this step is now deterministic, because everyone from the past time\n # step recovers\n new_s <- S - new_i\n new_r <- R + I\n\n cbind(new_s, new_i, new_r)\n}\n```\n:::\n\n\nThen we can run the simulation over a set of discrete time steps, T.\nNote again that T is now the number of infectious generation times, since time is rescaled relative to the models above.\n\n\n::: {.cell}\n\n```{.r .cell-code}\n####################################################################\n# Run over many steps\t\t\t\t\t\t\t\t\t\t\t #\n####################################################################\n\nmax_time <- 20\n# note here that the time step is one infectious generation time,\n# so 7 days from gamma above\nsims <- 1000\n\ns_mat <- matrix(S, 1, sims) # storage item for S for all simulations\ni_mat <- matrix(I, 1, sims) # storage item for I for all simulations\nr_mat <- matrix(R, 1, sims) # storage item for R for all simulations\nn_mat <- S + I + R\n\nfor (time_step in 2:max_time) {\n out <- chain_binomial_step(\n sims,\n s_mat[time_step - 1, ],\n i_mat[time_step - 1, ],\n r_mat[time_step - 1, ],\n beta\n )\n # update state\n s_mat <- rbind(s_mat, out[, 1])\n # update state\n i_mat <- rbind(i_mat, out[, 2])\n # update state\n r_mat <- rbind(r_mat, out[, 3])\n # update state -- note population size isn't changing, but this could be\n # updated with births/deaths\n n_mat <- rbind(n_mat, out[, 1] + out[, 2] + out[, 3])\n}\n```\n:::\n\n\nWe can plot, but note again, that this is plotted in terms of epidemic generations on the X-axis, rather than days (as in the previous sections).\n\n\n::: {.cell .column-body}\n\n```{.r .cell-code}\n# put output in a data frame\nbin_df <- tibble(\n S = array(s_mat),\n I = array(i_mat),\n R = array(r_mat),\n N = array(matrix(n_mat)),\n time = rep(1:max_time, sims),\n iteration = rep(1:sims, each = max_time)\n)\n\nbin_df <- bin_df %>%\n pivot_longer(\n # Don't pivot the time column\n cols = -c(time, iteration),\n names_to = \"state\",\n values_to = \"number\"\n ) %>%\n mutate(state = factor(state, levels = c(\"S\", \"I\", \"R\", \"N\")))\n\n# fix this color\nsir_colors <- c(S = \"#1f77b4\", I = \"#ff7f0e\", R = \"#FF3851\")\n\nbin_df %>%\n mutate(iteration = as.factor(iteration)) %>%\n filter(state %in% c(\"S\", \"I\", \"R\")) %>%\n ggplot(aes(\n x = time,\n y = number,\n group = interaction(iteration, state),\n color = state\n )) +\n geom_line(linewidth = 1.5, alpha = .1) +\n scale_color_manual(values = sir_colors) +\n labs(\n x = \"Time\",\n y = \"Number\",\n color = \"State\"\n ) +\n guides(color = guide_legend(override.aes = list(alpha = 1))) +\n theme(legend.position = \"top\")\n```\n\n::: {.cell-output-display}\n![](r-session-04_files/figure-html/unnamed-chunk-23-1.png){width=100%}\n:::\n:::\n\n", "supporting": [], "filters": [ "rmarkdown/pagebreak.lua"