In this practical we will examine the effect of changing the time step on the model outputs. Remember that in the lecture we said that we increase the accuracy of simulations when we use small time steps, and we reduce the accuracy (but speed up the simulation) when we use large time steps.
So far, we have written code to model difference equations, where events happen at fixed intervals. For some biological populations such as insects reproducing on an annual cycle this framework is very appropriate. However, in general we are interested in events that could happen at any point, and importantly once an event has happened it can affect the rate at which new events happen – for instance, when the first infected animal in a large population succeeds in infecting a second, the rate of infection approximately doubles. If we use large time steps we may miss this effect of changing transmission rates between one time step and the next, whereas if we use very small timesteps, we will be accurate but many steps will contain no events, thus wasting computational time. In this practical we will see that the size of the time step affects our model.
To change the model to handle timesteps, we have to consider that by default, in the previous programs we were using a timestep of 1. In this case – of E. coli transmission in cattle from Practical 2-1 – these are measured in weeks, which means that if our rates were given to us in units of weeks, then we are stepping forward one week at a time. If our rates were given to us per day rather than per week we would be stepping forward day by day. If a rate is given to us in weeks, we could scale it down accordingly to get the rate per day. Or, we could scale it up to get the rate per month. So to change our program to step forward on a different time scale, we scale the parameters. Our parameters (the transmission and recovery rates) will help us determine the number of transmissions and recoveries that occur in one unit of time. To simulate stepping forward in larger of smaller amounts of time we scale the parameters by the new time step.
The task now is to edit 0201-run-SIS.R so that it is easy to
alter the timestep of the simulation – you can reuse
step_deterministic_SIS(). All you need to do is to write a
new script file which creates a new parameter, timestep, so
that you can scale transmission.rate and
recovery.rate to ecoli.transmission * timestep
and ecoli.recovery * timestep before passing them to
step_deterministic_SIS(). You also need to remember that
the vector of timesteps to loop through is now also wrong.
Make sure that you adjust it correctly so that it will give the correct
times for the plot. You can check it has by looking at
population.df$time after you’ve added it in to the data
frame and make sure it looks right to you. We can call this new version
0202-test-SIS.R.
timestep does not need to be a whole number, and it can
be less than 1. Start running your programming with a
timestep of 1, and then try reducing it to 0.1 and then
compare your outputs. If they appear unchanged it means the first
timestep was small enough to produce an accurate output. Remember you
can plot multiple graphs on top of each other to make it easier to
directly compare timesteps.
Once you’ve satisfied yourself that the simulation is currently
accurate, find out by how much you can increase the timestep safely
before the simulation becomes inaccurate. Running this code with a
timestep of 3 or 5 means that we are stepping forward in
chunks of time that are 3 or 5 times larger than our default time step.
Note that if you don’t use the parameters from E. coli
transmission in cattle, the values of the timestep at which is all goes
wrong may be different, so we recommend sticking with those values.
Demonstrate the instabilities at high timesteps in a report, with some text explaining when things break down and the simulation is no longer reliable, and when, conversely, there is no advantage to reducing the timesteps further.
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.
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.
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 as well if they are there).
Remember that interacting like this through GitHub to help each other will count as your engagement marks for the course.