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. The other vignettes and the odin & monty book 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. \(\beta\) is the infection rate, \(\gamma\) 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