Code your first stochastic model. (2024)

  • 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 package adaptivetau
  • 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

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

## 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

Code your first stochastic model. (1)

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:

## $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

Code your first stochastic model. (2024)

References

Top Articles
Latest Posts
Recommended Articles
Article information

Author: Rev. Porsche Oberbrunner

Last Updated:

Views: 6145

Rating: 4.2 / 5 (53 voted)

Reviews: 92% of readers found this page helpful

Author information

Name: Rev. Porsche Oberbrunner

Birthday: 1994-06-25

Address: Suite 153 582 Lubowitz Walks, Port Alfredoborough, IN 72879-2838

Phone: +128413562823324

Job: IT Strategist

Hobby: Video gaming, Basketball, Web surfing, Book restoration, Jogging, Shooting, Fishing

Introduction: My name is Rev. Porsche Oberbrunner, I am a zany, graceful, talented, witty, determined, shiny, enchanting person who loves writing and wants to share my knowledge and understanding with you.