odin2 implements a high-level language for describing and implementing ordinary differential equations and difference equations in R. It provides a “domain specific language” (DSL) which looks like R but is compiled directly to C++, using dust2 to solve your system and to provide an interface to particle filters. You can then use monty to fit your models using MCMC.

This vignette jumps through a few of the core features of odin2 and ways that you might use it with dust2 and monty. Other vignettes (when written!) will expand on topics covered here in more detail.

Discrete time stochastic SIR model

A simple definition of the SIR model is:

$$\begin{align*} \frac{dS}{dt} &= -\beta \frac{SI}{N} \\ \frac{dI}{dt} &= \beta \frac{SI}{N} - \gamma I \\ \frac{dR}{dt} &= \gamma I \\ \end{align*}$$

where S is the number of susceptibles, I is the number of infected and R is the number recovered; the total population size N = S + I + R is constant. β is the infection rate, γ is the recovery rate.

Discretising this model in time steps of width dt gives the following update equations for each time step:

$$\begin{align*} S_{t+1} &= S_t - n_{SI} \\ I_{t+1} &= I_t + n_{SI} - n_{IR} \\ R_{t+1} &= R_t + n_{IR} \end{align*}$$

where

$$\begin{align*} n_{SI} &\sim \mathrm{Binomial}(S, 1 - e^{-\beta \frac{I}{N} \cdot dt}) \\ n_{IR} &\sim \mathrm{Binomial}(I, 1 - e^{-\gamma \cdot dt}) \end{align*}$$

Here is this system, as a stochastic compartmental model:

gen <- odin2::odin({
  p_IR <- 1 - exp(-gamma * dt)
  N <- parameter(1000)

  p_SI <- 1 - exp(-(beta * I / N * dt))
  n_SI <- Binomial(S, p_SI)
  n_IR <- Binomial(I, p_IR)

  update(S) <- S - n_SI
  update(I) <- I + n_SI - n_IR
  update(R) <- R + n_IR

  initial(S) <- N - I0
  initial(I) <- I0
  initial(R) <- 0

  beta <- parameter(0.2)
  gamma <- parameter(0.1)
  I0 <- parameter(10)
})

This step generates C++ code for the model and compiles it; it will take a few seconds.

Once we have our system, we can pass it to [dust2::dust_system_create] to create and start simulating it. Our system above has defaults for its parameters (N, beta, gamma, and I0) so we can initialise with almost no arguments:

pars <- list(beta = 0.2, gamma = 0.1, I0 = 10, N = 1000)
sys <- dust2::dust_system_create(gen(), pars, n_particles = 10)

By default the system will start at time 0 and with dt = 1. We can simulate 10 random epidemics starting from our initial conditions:

dust2::dust_system_set_state_initial(sys)
time <- 0:100
y <- dust2::dust_system_simulate(sys, time)
matplot(time, t(y[2, , ]), col = "#00000055", lty = 1, type = "l",
        xlab = "Time", ylab = "Number of infecteds")

as this system is stochastic, each trajectory will be different.

- Discrete time stochastic SIR model

Getting started with odin2

Discrete time stochastic SIR model