knitr::opts_chunk$set( collapse = TRUE, comment = "#>" )
library(epidmod)
epidmod is a system to simulate three typical stochastic epidemic models: stochastic SIR, SEIR, SEIQHRF model which provides a great deal of control and flexibility in the model itself and arrangement of information output. epidmod does not show the numerical calculations and What it does provide is a basic idea for stochastic simulation.The aim of this document is to provide an introduction to the usage of fundamental functions underlying the epidmod package: rSIR, rSEIR, rSEIQHRF, the corresponding print & plot functions to show the simulation outputs and other S3 methods to extract and study information from the simulation. What can be achieved and how to achieve it using epidmod are also demonstrated by several examples below.
Usage of rSIR and rSEIR
rSIR and rSEIR are designed to do the simulation of stochastic SIR and SEIR model with matrix outputs. The SIR model is the most basic one and SEIR, SEIQHRF models are all its extensions. When applying the functions, users must be very careful about the parameters. The below codes demonstrate the correct usage with random selected parameters.
model1 <- rSIR (N0 = 1000, I0 = 50, S0 = 950, days = 300, pars = c(1/10, 1, 1/5, 0.1, 0.1)) model2 <- rSEIR(N0 = 1000, I0 = 50, S0 = 949, R0 = 0, E0 = 1 , days = 300, pars = c(1/10, 1, 1/4, 1/5, 0, 0.1, 0.1))
The objects returned by the rSIR and rSEIR function are two complicated matrix to store all the values in each groups during the simulation and S3 generic methods are used for presenting better view of outputs (will be mentioned in next section).The first several arguments represent the set initial value of number of people in different groups and the total simulation time (the unit to measure time in this system is day).It is obvious that this two function are very similar and they both have an argument pars to store the needed parameters of transition rates. Now we introduce the meaning of each parameters in par:
Parameters stored in pars of rSIR
- a.rate rate of arrival of people into the population
- i.rate individual effective contact rate
- r.rate recovery time for each infected individual is 1/ r.rate
- pI probability that an arrival is infected
- pIm probability that an infected person is immune after recovery
Parameters stored in pars of rSEIR
- a.rate rate of arrival of people into the population
- e.rate individual effective contact rate
- i.rate the incubation period is 1/i.rate days
- r.rate recovery time for each infected individual is 1/ r.rate
- pE probability that an arrival is exposed.
- pSR probability that an exposed individual self-recover and turns susceptible.
- pIm probability that an infected person is immune after recovery
Most of the parameters have the same meaning, however, the user should be careful about the i.rate when applying the function. Due to the different assumptions in these two models, i.rate is the object to control the incubation time.
Usage of rSEIQR and rSEIQHRF
In this package, SEIQHRF is the most complicated model and it contains seven groups and much more parameters. The simulation is harder than previous ones but applying the function rSEIQHRF is still very easy to handle. The usage of rSEIQR is basically the simplified version of rSEIQHRF. The examples are below:
model3 <- rSEIQR(N0 = 100, S0 = 99, E0 = 1 ,I0 = 0, Q0 = 0, R0 = 0, days = 100, pars = c(1/12, 1, 0.5, 0.15, 0.04, 0, 0, 1)) para <- c( 1/10, 10, 0.01, 0.02, 0.2, 0.01, 1/20, 1/15, 1/30, 1/1000, 0.01, 0.2, 0.01, 0.02, 0.03) model4 <-rSEIQHRF(N0 = 500, S0 = 497, E0 = 0, I0 = 3, Q0 = 0, H0 = 0, R0 = 0, F0 = 0, days = 365, pars = para )
Like previous functions, arguments are divided in two parts, the first half is still the initial value with more compartments and setting time, the other is about the stored parameters pars. Users should be really careful about the input parameters. The representations are below:
Parameters stored in pars of rSEIQR
- a.rate rate of arrival of people into the population
- e.rate individual effective contact rate
- i.rate the incubation period is 1/i.rate days
- q.rate infected individual self-quarantining rate at time t
- r.rate recovery time for each infected individual is 1/ r.rate
- pE probability that an arrival is exposed.
- pSR probability that an exposed individual self-recover and turns susceptible.
- pIm probability that an infected person is immune after recovery
Parameters stored in pars of rSEIQHRF
- a.rate rate of arrival of people into the population
- e.rate individual effective contact rate
- i.rate the incubation period is 1/i.rate days
- h1.rate infected individual requiring hospitalization rate at time t
- q.rate infected individual self-quarantining rate at time t
- h2.rate self-isolated individual requiring hospitalization rate at time t
- r1.rate recovery time for each infected individual is 1/ r1.rate
- r2.rate recovery time for each individual who accepts hospitalization is 1/ r2.rate
- r3.rate recovery time for each individual who is self-quarantined is is 1/ r3.rate
- f.rate individual who accepts hospitalization case fatality rate at time t
- pE probability that an arrival is exposed
- pSR probability that an exposed individual self-recover and turns susceptible
- pIim probability that an infected person is immune after recovery
- pQim probability that a quarantined person is immune after recovery
- pHim probability that an individual who requires hospitalization recovers and turns susceptible
Print & Plot method
The returns in above functions are all huge matrix with stored value and that's inappropriate for presenting. To solve this problem, I introduce the S3 print and plot method to help present the outputs.Let use the previous model1 as the example to show the output:
print(model1)
The print method firstly introduces the type of the model and then a simulation summary would be presented which includes the maximum infection rate during the whole simulation and the needed time to reach that maximum. The motivation about this summary is to calculate out the worst situation with given inputs and remaining time for people to deal with the epidemic.Then comes to the model parameters part, where all parameters involved in the simulation would be presented for calibration of inputs (this would be especially useful when applying the rSEIQHRF function). Last part is to show the number of people in each groups at the end of simulation and this would give users a very intuitive feeling about what is simulated at the end.
plot(model1)
The plot method produces a basic time series plot with coordinate time and number of people. Different groups of curves are marked with different colors, the line type and point shape of the curve with legends on top right (In order to ensure that the image is easy to distinguish, in the SEIQHRF model, only color but line types look the same). Basically, this is a very fundamental plot method but still provide the users a very clear look about the trend of the entire simulation and help them to understand the meaning of each parameter. For example, if we make some adjustment on the parameter as below, the plot will be like:
model1 <- rSIR (N0 = 1000, I0 = 50, S0 = 950, days = 300, pars = c(1/10, 1, 1/5, 0.1, 0.1)) model5 <- rSIR (N0 = 1000, I0 = 50, S0 = 950, days = 300, pars = c(1/10, 1.5, 0.1, 0.1, 0.1)) plot(model5)
Comparing Figure11 with Figure12, we can clearly see that the peak of infection is much higher than before and the number of susceptible goes to a lower level. This is reasonable since we select a higher contact rate $(\beta = 1\rightarrow\beta = 1.5)$ which indicates the disease in model5 is more contagious. Figure12 also suggests it require more time for recovered people to reach stable level e which is due to the lower recovery rate $(\delta = 0.2 \rightarrow \delta = 0.1)$ in model5 and this indicates the average recovery time increases from $ 1/0.2 =5 $ to $ 1/0.1 = 10$. If we keep the parameter fixed and change the input number of people in each groups, we will find out more about the impact of initial value in each stochastic model.
plot(model2)
The parameter settings in model1 and model2 are basically the same except that there is a parameter $\beta = 0.25$ (in SEIR model) representing $1/0.25 = 4$ incubation period in model2. Comparing Figure 13 with Figure 11, it illustrates that the infected group in Figure 11 is divided into two new groups: Exposed and Infected in Figure 13 which is coincident with the meaning of incubation period(explained in Section2.3). The new phenomenon in Figure 13 is that there seems less recovered people and more susceptible people in the later stages of the simulation for SEIR model. This could be treated as the effect of adding a new compartment on the original model or explained by the interaction effect between the Infected group and Exposed group.
Other S3 methods
In this package, there are other functions to extract information from the simulated results and they will be demonstrated one by one.
1.IR_info
This function produces some basic information similar with that in print but with a specific infected proportion I(t)/N(t) and it is useful when users would like find out how is the simulation with a given specific proportion of infected people. The original intention of designing it is to make the whole simulation more flexible and it also show the needed time to reach that given infection rate and the number of people in each group at that time. User should be careful about the input infection should be equal or lower than the maximum infected proportion (could be found in output of print).
IR_info(model1, 0.3)
2. hist
This function is built to study the distribution of the number of people in each groups at a given time and total number of simulation. The core idea is to simulate the same model with fixed input and parameters many times then choose a selected input time and collect data about the number of people in each group at that time in each simulation to produce the histogram for each group. This is very helpful since if we can find some pattern from the distribution of groups with more concerns(for example, the infected in all models or the group of people who requires hospitalization), it will help people better prepared for restricting and controlling the disease. The second argument is the input time(should be less than the total simulation time) and third argument is to determine the total number of simulation. The output of this function is multiple histograms.
model5 <- rSIR (N0 = 1000, I0 = 200, S0 = 800, days = 300, pars = c(1/10, 1, 1/5, 0.1, 0.1)) hist( model5, 70, 100)
3. afterpeak_info
This function returns the peak information about the simulation and after peak, the time until the infected reduce to 50%, 30% and 10% of the maximum infected. It will help the users to study the after peak information, estimate time when the epidemic will come to an end and better assess the impact of the epidemic. The usage of this function is quite easy with only one input but it is very effective when predicting how long the disease will affect society negatively and estimating the revival time of society and economics.
afterpeak_info(model1)
