In IBAHCM/RPiR: Reproducible Programming in R

library(learnr)

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

Overview

In this practical we separate out the model from the simulation functionality more cleanly. In the process we make the previous set of code safer to use (and more elegant!).

Background

When we are programming we often make quick additions to code, as we did in the last practical but it is good practice, if we want to keep them, to make sure they are included in the code in a way that makes it clear what's going on. We had to know how the SIS model works to know that we could get away with multiplying the parameters with timestep. This use of a low-level detail of how the model runs might be risky -- we may later change the model and not realise this hack no longer works.

Tasks

In the previous practical we multiplied the transmission and recovery rates up (or down) to reflect changes in timestep. However, transmission.rate and recovery.rate were not really changing, and we want the code to reflect that rather than abusing it to do something it wasn't intended to do. Instead add a new argument, timestep, into step_deterministic_SIS() -- we could call it timestep_deterministic_SIS(). Rather than submit the product transmission.rate * timestep as an argument to the simulation function, we can now separately pass the real rate parameters and the timestep separately. This rescaling is now done within the function timestep_deterministic_SIS() in 0203-deterministic-SIS.R. We might, for instance, create two new variables, perhaps called effective.transmission.rate and effective.recovery.rate (or something shorter!) inside that function. These new variables are just our original rates rescaled to match the new time scale.

Now that we have a timestep in the step function, it makes sense to have the time in there too, and that can be updated inside the function. Indeed, the data frame itself that is passed back and forth can have an extra column. This should start in the new report script that you are writing (perhaps 0203-run-SIS.R), where it will be initialised with all three columns (time, susceptibles and infecteds):

population.df <- data.frame(time = start.time,
                            susceptibles = initial.susceptibles,
                            infecteds = initial.infecteds)

And the timestep_deterministic_SIS() function can update the time along with the number of susceptibles and infecteds. Here, it makes the code slightly more elegant but it might be critical if, for instance, we wanted to adjust the timestep as we went along, which we will need to do later. Now we're not making the vector of timesteps in advance, we can use a while() loop to decide when we're finished. For instance:

latest.population <- population.df
while (latest.population$time < end.time) {
  latest.population <- timestep_deterministic_SIS(latest = latest.population,
                                                  transmission.rate = ecoli.transmission, 
                                                  recovery.rate = ecoli.recovery, 
                                                  timestep = this.timestep)
  population.df <- rbind(population.df, latest.population) 
}

latest.population$time is the time that we are currently at, which we extract from the data frame after each timestep.

Running the code

Run this code and check that it gives the same answers as in the previous practical. Don't forget to check that your function is using only arguments passed to it by running findGlobals().

Report

Demonstrate that the code is the same in a report by running both step_deterministic_SIS() and timestep_deterministic_SIS() with a couple of different values of timestep (using source() to load in both functions from their respective files).

GitHub

As with the other exercises, we want you to get a couple of other people to check your code and make sure it works for them, and we want you to check other people's code too.

Getting help with your code

Once you're happy with your code, commit your changes using the Git pane in RStudio (just commit the R files - not, for instance, any html files that were created by generating reports). Don't forget then to push the changes to GitHub and check on the website that it contains your new code. Notify your partners that you have something for them to check. If you're not sitting next to them, then you can create an issue in your repository asking for their review - you can contact them by tagging them in the message with @theirusername, and they should receive an email and a notification on the website. They should then create an issue, telling you if they had any problems running your version of the practical.

Checking other people's code

Likewise, when you hear from someone else that a repo is available to check, then if it's a project you have already downloaded, you should open that project in RStudio, go to the Git pane, and click on the Pull button to get the updates they have pushed to GitHub (there should be a subtle message saying that updates are available near the button). If you're checking a new person's work, or a new project, refer to the instructions in Practical 2-1. Once you've checked that the new code has downloaded, you can then run their code and make sure it works. Open an issue on GitHub to say if it works, or explaining what went wrong if it doesn't (it might help to also do this in person if they are there).

Interacting like this through GitHub to help each other will count as your engagement marks for the course.

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