In IBAHCM/RPiR: Reproducible Programming in R

library(learnr)

knitr::opts_chunk$set(error = TRUE)
set.seed(123)

Overview

In this practical we develop and run our first stochastic version of the SIS model.

Background

We saw in the previous practicals the impact of including randomness into the exponential model with births and deaths. Having developed the programming framework, it is now time to investigate the impact of randomness in the simplest epidemiological model -- the SIS model.

Tasks

In this practical we will write a stochastic version of the SIS model. You have been provided with run_simulation() in the RPiR package, which allows some extra functionality over run_simple() to run simulations using the older deterministic models. Pay attention in particular to writing the function, where the most important changes need to be made.

1. Write a function

First, you need to write a new function, timestep_stochastic_SIS(), adapted from timestep_stochastic_birth_death() to handle the SIS model rather than exponential growth.

The function should be defined like this:

timestep_stochastic_SIS <- function(latest, 
                                    transmission.rate, 
                                    recovery.rate, 
                                    timestep) {
  # Some code

}

It proceeds much as before with the deterministic SIS model, calculating the effective transmission rate (beta) and recovery rate (sigma) at this timestep. Then as with the stochastic birth death model, you need to work out whether the individual rate of infection or recovery per timestep is greater than 1, in which case it cannot be used as a probability. In this case, the actual individual infection rate is effective.transmission.rate $\times \frac{I}{N}$ and the individual recovery rate is effective.recovery.rate.

Now for the key step... you need to sample from the binomial distribution (as described in the lecture) to determine the number of new infecteds and susceptibles using rbinom() with the probability of "success" set to the individual rates detailed above. Then update the populations, and if there are no infecteds left, then the epidemic is over, so set is.finished to TRUE.

is.finished <- next.infecteds == 0

Then append the new data to the data frame as usual, and return a list with the updated population data frame and a variable that tells us whether we are finished or not.

list(updated.pop = next.population, end.experiment = is.finished)

2. Write a demo

Now write a new demo called d0304_run_SIS, to demonstrate the SIS model. It might be useful to adapt the code you wrote in d0303_run_birth_death.

3. Check your function

What happens when you do this when the time step is 1? You should get an error message saying the timestep is too large. You may still get a plot though! This is because R is plotting the previous value of populations held in its memory. This is why it's important to be aware of your workspace and clear your previous history when you're working on a new project (don't use rm(list=ls()) though!).

Remember that generating a report starts a new copy of R and so is a much better check of whether the script can run on its own.

Clear what is held in R's workspace in RStudio. Now rerun with $R_0 = 15$ and confirm that you don't get an output when the error message shows. Remember that clearing your workspace and trying again is often a very useful trick when you get unexpected results.

Report

1. Compare stochastic and deterministic models

When you are happy with the structure of the code, run it a few times with a transmission rate of 0.3, a recovery rate of 0.1, and 2 initial infecteds from a population size of 100. Demonstrate the variability in the results of the stochastic model.

You should also overlay the deterministic results from Practical 2-3. Noting that run_simulation() will work with this as well since it automatically detects that it just returns a dataframe, and behaves appropriately:

final.populations <- run_simulation(timestep_deterministic_SIS, 
                                    initial.populations,
                                    end.time,
                                    transmission.rate = ecoli.trans,
                                    recovery.rate = ecoli.recov, 
                                    timestep = working.timestep)

2. Vary the rate of transmission and the initial number of infecteds

Try varying the transmission rate (beta) -- this will change $R_0$ -- and the initial number of infecteds and see how the likelihood of early extinction changes.

3. Vary the total population size

Try changing the total population size and determine the impact on the observed variability.

Check it works

As with previous exercises, you need to check that everything works correctly -- that the package installs, and the demos and help files work and you can generate reports from the demos -- and then we want you to get a couple of other people in your subgroup to check your code and make sure it works for them, and we want you to check other people's code too. We describe how to do this for packages in GitHub in Practical 3-1 (also under Check it works) if you're uncertain.

Remember, interacting like this through GitHub to help each other will count as most of your engagement marks for the course.

IBAHCM/RPiR documentation built on Jan. 12, 2023, 7:41 p.m.