Background

Preliminaries

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

  • draw samples from the distributions of the final size and the duration of the epidemic (see below for definitions of these quantities);
  • be able to simulate from an epidemic model where the infectious period follows a Weibull distribution;
  • be able to simulate from a stochastic Susceptible-Exposed-Infective-Removed (SEIR) epidemic model.

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.

Definitions

In this section we give definitions of two quantities for which samples from their distribution will be drawn.

The threshold parameter \(R_0\)

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

The final size of an epidemic is the total number of initial susceptibles who contracted the disease by the end of the outbreak.

The duration of the epidemic

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.

Exercises

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")

Exercise 1

Go through the lines of each function and make sure that you follow the logic behind.

Exercise 2

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)\).

Exercise 3

Modify the existing functions in simulation.R to record the final size and the duration of the epidemic as part of the functions’ output.

Exercise 4

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:

  1. Simulate a realisation from the epidemic model;
  2. 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.

Exercise 5

Repeat the above exercise but derive, by simulation, the distribution of the duration of the epidemic instead of the final size.

Exercise 6

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.

Exercise 7

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.