In the first lecture we introduced the concept of simulation by which we meant the procedure of producing a realisation of the model, or in other words a possible outcome. In this lab session we will be modifying existing R functions and by the end of this session, among other things, we will
Start by downloading the file simulation.R
from here
and save it in your workspace. This file simulation.R
contains four different functions which will enable us to simulate
realisations from various stochastic epidemic models:
simSIR.Markov
: This function produces realisations
from a Markovian SIR model with infection and removal rate, \(\beta/N\) and \(\gamma\) respectively. The procedure (as
described in Section 2 in the lecture) is the following: i) first
simulate the time to the next event and then ii) decide the
type of the event (infection or removal).
simSIR.Markov.alternative
: This function also
produces realisations from a Markovian SIR model with infection and
removal rate, \(\beta/N\) and \(\gamma\) respectively but with a slightly
different algorithm to the one described above. The procedure here is as
follows: i) generate a possible time to the next infection and the
possible time to the next removal; then ii) the event which happens
first determines the type of the next event.
simSIR.Non.Markov.constant
and
simSIR.Non.Markov.gamma
: These two functions allows us to
simulate realisations from non-Markov SIR models; with constant and
Gamma infectious period respectively. The procedure is similar to the
one used in simSIR.Markov.alternative
based on Section 5 in
Lecture 1.
In this section we give definitions of two quantities for which samples from their distribution will be drawn.
An appealing feature of an epidemic model is that it embodies a threshold parameter which can be utilised as a severity measure of the outbreak. Stochastic models such as epidemics and branching processes typically generate bimodal realisations where an epidemic may or may not die out quickly, depending on the value of a threshold parameter \(R_0\) often referred to as the basic reproduction number.
This is the most important parameter in epidemic theory and it is defined as the expected number of infections generated by a typical infective in an infinite susceptible population.
In the case of a homogeneously mixing stochastic epidemic, it holds that \[R_0 = \beta\mathbb{E}[I]\] where \(\mathbb{E}[I]\) denotes the expected infectious period. For instance, in the case of the Markovian SIR model, \[R_0 = \frac{\beta}{\gamma},\] since the individuals remain infectious for an average period of length \(1/\gamma\).
The final size of an epidemic is the total number of initial susceptibles who contracted the disease by the end of the outbreak.
Another quantity of interest is the duration of the epidemic and this is defined as the time that elapsed from the initial infection until the time of the last removal.
Download the file simulation.R
from the module’s
website and load all the functions into R by using the command
source
:
source("simulation.R")
Go through the lines of each function and make sure that you follow the logic behind.
Simulate realisations from a Markovian SIR using the function
simSIR.Markov
and make sure that you understand the output.
You may assume that size of the population size is \(N=21\) (i.e. 20 susceptibles and 1 initial
infective). In addition, you could try different values for \((\beta,\gamma)\), e.g. \((0.9, 1), (2, 1)\) and \((4, 1)\).
Modify the existing functions in simulation.R
to record
the final size and the duration of the epidemic as
part of the functions’ output.
Derive a simulation-based estimate of the distribution of the final size of a Markovian SIR model for different values of \(R_0\), e.g. \(R_0 = 0.9\), \(R_0 = 1.5\) and \(R_0 = 4\). Furthermore, do the same for the non-Markovian models, e.g. for a constant and a Gamma infectious period. Hint: You may find it useful to write a loop which will iterate the following steps for a number of times:
- Simulate a realisation from the epidemic model;
- Store the final size
At the end you should have a collection of final sizes for which then you can plot a histogram as your estimate of the true distribution of the final size.
Repeat the above exercise but derive, by simulation, the distribution of the duration of the epidemic instead of the final size.
Write a function in R to simulate from a non-Markovian stochastic epidemic model where the infectious period is assumed to follow a Weibull distribution. Hint: The probability density function (pdf) of the Weibull distribution is as follows: \[f(x) = (a/b) (x/b)^{(a-1)} \exp(- (x/b)^a), \qquad x > 0, \quad a > 0, \quad b > 0\] Type
?Weibull
to find out how to simulate from a Weibull distribution.
Write a function to simulate from an epidemic model which involves a fixed latent period, i.e. write a function to simulate from a stochastic SEIR model.