- Objectives
- Design your favourite model
- Write the transition table
- Implement the model in R
- Going further
- Observation time
- Computing incidence
- Design your favourite model with pen & paper.
- Write down a table with all transitions and jump intensities for this model.
- Implement the model in
R
using the packageadaptivetau
- Simulate and plot stochastic trajectories to investigate on the dynamics of your model.
Using only pen & paper, design a compartmental model with boxes and arrows. You can design either:
- A model you would like to work with for a particular disease in your research. In this case, try to choose one not too complicated, say with less than 10 compartments, otherwise it will be tricky to debug.
- The SEITL model, which is one of the 5 models used to analyse the two-wave flu outbreak in Tristan da Cunha. You can use the information provided in this page to design the SEITL model.
Include labels for the rates of flow between compartments and define all the symbols and parameters you use. You should also give a brief description of what the model does as well as some guess values for the parameters and initial conditions.
Once you have designed your model, write down the associated transition table. The table should contain 3 columns:
- Event: describe the type of transition (e.g.infection)
- Transition: change in the state vector corresponding to the transition (e.g. \((s, i, r)\to(s-1, i+1, r)\))
- Jump intensity: intensity at which the transition occur (e.g. \(\beta si/N\))
If you chose the SEITL model you check your table with our solution.
The adaptivetau package can be installed with
Once installed, it is loaded with
The adaptivetau
package uses a different syntax from the deSolve
package. Instead of providing a function to calculate the rates of change at each time point, one specifies a list of transitions and their rates (jump intensities). Examples for how this is done can be found in the adaptivetau vignette.
For the SIR model, we could write
SIR_transitions <- list( c(S = -1, I = 1), # infection c(I = -1, R = 1) # recovery )SIR_rateFunc <- function(x, parameters, t) { # define model parameters in term of the natural parameters beta <- parameters["R0"]/parameters["D_inf"] nu <- 1/parameters["D_inf"] # create temporary variables for states to simplify the writing of the rates below S <- x["S"] I <- x["I"] R <- x["R"] N <- S + I + R # return rates return(c( beta * S * I / N, # infection nu * I # recovery )) }
You can now re-use and modify the code of the SIR model to implement your model in R
. If you chose the SEITL model and you need more help, you can use and complete this partial solution.
To run the stochastic model, we then use the ssa.adaptivetau
function, which takes a vector of initial conditions, the list of transitions and rate function, a named vector of parameters, and the final time (with simulations starting at time 0).
parameters <- c(R0 = 5, D_inf = 1)run <- ssa.adaptivetau(init.values = c(S = 999, I = 1, R = 0), transitions = SIR_transitions, rateFunc = SIR_rateFunc, params = parameters, tf = 10)head(run)
## time S I R## [1,] 0.0000000 999 1 0## [2,] 0.4173075 998 2 0## [3,] 0.4178569 997 3 0## [4,] 0.4363569 996 4 0## [5,] 0.4525435 995 5 0## [6,] 0.4801278 995 4 1
The simplest way to plot the trajectories is using plot
. To plot the output of the stochastic SIR run above, we first convert it to a data frame (ssa.adaptivetau
returns a matrix) using data.frame
We can then plot the number of infected against time using
Spend some time playing with your stochastic model by changing the parameter values and see how it affects the epidemic. In particular, try a value of \(R_0\) just above 1 to see the early stochastic extinction.
Observation time
Unlike the function ode
from the deSolve
package, ssa.adaptivetau
does not produce output at specific times, but every time an event happens. However, in reality, we observe epidemics at specific points in time, for instance once per day.
To get the output at chosen times, we can use approx
:
# get output at days 1, 2, ... , 10run_I_times <- approx(x = run_df$time, y = run_df$I, xout = 1:10, method = "constant")print(run_I_times)
## $x## [1] 1 2 3 4 5 6 7 8 9 10## ## $y## [1] 16 325 363 170 77 32 16 7 3 1
Write a function that will apply the approx
function to all the variables returned by ssa.adaptivetau
and construct a data frame with model output at the desired times.
Computing incidence
You may have noticed that the Tristan da Cunha epidemic dataset corresponds to daily incidence (number of new cases infected per day) rather than prevalence (number of cases infected over time). This is a common feature of epidemiological data: it is much easier to report new cases as they present to their doctor than counting the number of cases currently infected on a given day.
The compartemental model gives prevalence rather than incidence. To obtain the daily incidence from our model, we need to create a \(6^\mathrm{th}\) state variable - called \(\mathrm{Inc}\) - to track the daily number of new cases. Assuming that new cases are reported when they become symptomatic and infectious, we need to modify the transition corresponding to the onset of infectiousness in the SEITL model:
\[(s, e, i, t, l, inc)\to(s, e-1, i+1, t, l, inc+1)\]
Modify your stochastic model to include this new transition and retrieve the daily incidence from the output of ssa.adaptivetau
. You will have to also make use of the function diff
, which computes the difference between two consecutive values of a vector
Top: Index Previous: R introduction Next: My first fitmodel