In its simplest form, a SIR model is typically written in continuous time as:
$$ \frac{dS}{dt} = - \beta \frac{S_t I_t}{N_t} $$
$$ \frac{dI}{dt} = \beta \frac{S_t I_t}{N_t} - \gamma I_t $$
$$ \frac{dR}{dt} = \gamma I_t $$
Where β is an infection rate and γ a removal rate, assuming ‘R’ stands for ‘recovered’, which can mean recovery or death.
For discrete time equivalent, we take a small time step t (typically a day), and write the changes of individuals in each compartment as:
$$ S_{t+1} = S_t - \beta \frac{S_t I_t}{N_t} $$
$$ I_{t+1} = I_t + \beta \frac{S_t I_t}{N_t} - \gamma I_t $$
Rt + 1 = Rt + γIt
The discrete model above remains deterministic: for given values of the rates β and γ, dynamics will be fixed. It is fairly straightforward to convert this discrete model into a stochastic one: one merely needs to uses appropriate probability distributions to model the transfer of individuals across compartments. There are at least 3 types of such distributions which will be useful to consider.
This distribution will be used to determine numbers of individuals leaving a given compartment. While we may be tempted to use a Poisson distribution with the rates specified in the equations above, this could lead to over-shooting, i.e. more individuals leaving a compartment than there actually are. To avoid infecting more people than there are susceptibles, we use a binomial distribution, with one draw for each individual in the compartment of interest. The workflow will be:
determine a per-capita probability of leaving the compartment, based on the original rates specified in the equations; if the rate at which each individual leaves a compartment is λ, then the corresponding probability of this individual leaving the compartment in one time unit is p = 1 − e−λ.
determine the number of individuals leaving the compartment using a Binomial distribution, with one draw per individual and a probability p
As an example, let us consider transition S → I in the SIR model. The overall rate at which this change happens is $\beta \frac{S_t I_t}{N_t}$. The corresponding per susceptible rate is $\beta \frac{I_t}{N_t}$. Therefore, the probability for an individual to move from S to I at time t is $p_{(S \rightarrow I), t} = 1 - e^{- \beta \frac{I_t}{N_t}}$.
Poisson distributions will be useful when individuals enter a compartment at a given rate, from ‘the outside’. This could be birth or migration (for S), or introduction of infections from an external reservoir (for I), etc. The essential distinction with the previous process is that individuals are not leaving an existing compartment.
This case is simple to handle: one just needs to draw new individuals entering the compartment from a Poisson distribution with the rate directly coming from the equations.
For instance, let us now consider a variant of the SIR model where new infectious cases are imported at a constant rate ϵ. The only change to the equation is for the infected compartment:
$$ I_{t+1} = I_t + \beta \frac{S_t I_t}{N_t} + \epsilon - \gamma I_t $$
where:
individuals move from S to I according to a Binomial distribution $\mathcal{B}(S_t, 1 - e^{- \beta \frac{I_t}{N_t}})$
new infected individuals are imported according to a Poisson distribution 𝒫(ϵ)
individual move from I to R according to a Binomial distribution ℬ(It, 1 − e−γ)
This distribution will be useful when individuals leaving a compartment are distributed over several compartments. The Multinomial distribution will be used to determine how many individuals end up in each compartment. Let us assume that individuals move from a compartment X to compartments A, B, and C, at rates λA, λB and λC. The workflow to handle these transitions will be:
determine the total number of individuals leaving X; this is done by summing the rates (λ = λA + λB + λC) to compute the per capita probability of leaving X (p(X → ...) = 1 − e−λ), and drawing the number of individuals leaving X (n(X → ...)) from a binomial distribution n(X → ...) ∼ B(X, p(X → ...))
compute relative probabilities of moving to the different compartments (using i as a placeholder for A, B, C): $p_i = \frac{\lambda_i}{\sum_i \lambda_i}$
determine the numbers of individuals moving to A, B and C using a Multinomial distribution: n(X → A, B, C) ∼ ℳ(n(X → ...), pA, pB, pC)
odin
We start by loading the odin
code for a discrete,
stochastic SIR model:
## Core equations for transitions between compartments:
update(S) <- S - beta * S * I / N
update(I) <- I + beta * S * I / N - gamma * I
update(R) <- R + gamma * I
## Total population size (odin will recompute this at each timestep:
## automatically)
N <- S + I + R
## Initial states:
initial(S) <- S_ini # will be user-defined
initial(I) <- I_ini # will be user-defined
initial(R) <- 0
## User defined parameters - default in parentheses:
S_ini <- user(1000)
I_ini <- user(1)
beta <- user(0.2)
gamma <- user(0.1)
As said in the previous vignette, remember this looks and parses like R code, but is not actually R code. Copy-pasting this in a R session will trigger errors.
We then use odin
to compile this model:
## <odin_model> object generator
## Public:
## initialize: function (..., user = list(...), use_dde = FALSE, unused_user_action = NULL)
## ir: function ()
## set_user: function (..., user = list(...), unused_user_action = NULL)
## initial: function (step)
## rhs: function (step, y)
## update: function (step, y)
## contents: function ()
## transform_variables: function (y)
## engine: function ()
## run: function (step, y = NULL, ..., use_names = TRUE)
## Private:
## ptr: NULL
## use_dde: NULL
## odin: NULL
## variable_order: NULL
## output_order: NULL
## n_out: NULL
## ynames: NULL
## interpolate_t: NULL
## cfuns: list
## dll: discrete.deterministic.sir5c1f7b62
## user: I_ini S_ini beta gamma
## registration: function ()
## set_initial: function (step, y, use_dde)
## update_metadata: function ()
## Parent env: <environment: namespace:discrete.deterministic.sir5c1f7b62>
## Locked objects: TRUE
## Locked class: FALSE
## Portable: TRUE
Note: this is the slow part (generation and then compilation of C code)! Which means for computer-intensive work, the number of times this is done should be minimized.
The returned object sir_generator
is an R6 generator that
can be used to create an instance of the model: generate an instance of
the model:
## <odin_model>
## Public:
## contents: function ()
## engine: function ()
## initial: function (step)
## initialize: function (..., user = list(...), use_dde = FALSE, unused_user_action = NULL)
## ir: function ()
## rhs: function (step, y)
## run: function (step, y = NULL, ..., use_names = TRUE)
## set_user: function (..., user = list(...), unused_user_action = NULL)
## transform_variables: function (y)
## update: function (step, y)
## Private:
## cfuns: list
## dll: discrete.deterministic.sir5c1f7b62
## interpolate_t: NULL
## n_out: 0
## odin: environment
## output_order: NULL
## ptr: externalptr
## registration: function ()
## set_initial: function (step, y, use_dde)
## update_metadata: function ()
## use_dde: FALSE
## user: I_ini S_ini beta gamma
## variable_order: list
## ynames: step S I R
x
is an ode_model
object which can be used
to generate dynamics of a discrete-time, deterministic SIR model. This
is achieved using the function x$run()
, providing time
steps as single argument, e.g.:
## step S I R
## [1,] 0 1000.0000 1.000000 0.0000000
## [2,] 1 999.8002 1.099800 0.1000000
## [3,] 2 999.5805 1.209517 0.2099800
## [4,] 3 999.3389 1.330125 0.3309317
## [5,] 4 999.0734 1.462696 0.4639442
## [6,] 5 998.7814 1.608403 0.6102138
## [7,] 6 998.4604 1.768530 0.7710541
## [8,] 7 998.1076 1.944486 0.9479071
## [9,] 8 997.7198 2.137811 1.1423557
## [10,] 9 997.2937 2.350191 1.3561368
## [11,] 10 996.8254 2.583469 1.5911558
x_res <- x$run(0:200)
par(mar = c(4.1, 5.1, 0.5, 0.5), las = 1)
matplot(x_res[, 1], x_res[, -1], xlab = "Time", ylab = "Number of individuals",
type = "l", col = sir_col, lty = 1)
legend("topright", lwd = 1, col = sir_col, legend = c("S", "I", "R"), bty = "n")
The stochastic equivalent of the previous model can be formulated in
odin
as follows:
## Core equations for transitions between compartments:
update(S) <- S - n_SI
update(I) <- I + n_SI - n_IR
update(R) <- R + n_IR
## Individual probabilities of transition:
p_SI <- 1 - exp(-beta * I / N) # S to I
p_IR <- 1 - exp(-gamma) # I to R
## Draws from binomial distributions for numbers changing between
## compartments:
n_SI <- rbinom(S, p_SI)
n_IR <- rbinom(I, p_IR)
## Total population size
N <- S + I + R
## Initial states:
initial(S) <- S_ini
initial(I) <- I_ini
initial(R) <- 0
## User defined parameters - default in parentheses:
S_ini <- user(1000)
I_ini <- user(1)
beta <- user(0.2)
gamma <- user(0.1)
We can use the same workflow as before to run the model, using 10
initial infected individuals (I_ini = 10
):
## <odin_model> object generator
## Public:
## initialize: function (..., user = list(...), use_dde = FALSE, unused_user_action = NULL)
## ir: function ()
## set_user: function (..., user = list(...), unused_user_action = NULL)
## initial: function (step)
## rhs: function (step, y)
## update: function (step, y)
## contents: function ()
## transform_variables: function (y)
## engine: function ()
## run: function (step, y = NULL, ..., use_names = TRUE)
## Private:
## ptr: NULL
## use_dde: NULL
## odin: NULL
## variable_order: NULL
## output_order: NULL
## n_out: NULL
## ynames: NULL
## interpolate_t: NULL
## cfuns: list
## dll: discrete.stochastic.sir206fe457
## user: I_ini S_ini beta gamma
## registration: function ()
## set_initial: function (step, y, use_dde)
## update_metadata: function ()
## Parent env: <environment: namespace:discrete.stochastic.sir206fe457>
## Locked objects: TRUE
## Locked class: FALSE
## Portable: TRUE
set.seed(1)
x_res <- x$run(0:100)
par(mar = c(4.1, 5.1, 0.5, 0.5), las = 1)
matplot(x_res[, 1], x_res[, -1], xlab = "Time", ylab = "Number of individuals",
type = "l", col = sir_col, lty = 1)
legend("topright", lwd = 1, col = sir_col, legend = c("S", "I", "R"), bty = "n")
This gives us a single stochastic realisation of the model, which is of limited interest. As an alternative, we can generate a large number of replicates using arrays for each compartment:
## Core equations for transitions between compartments:
update(S[]) <- S[i] - n_SI[i]
update(I[]) <- I[i] + n_SI[i] - n_IR[i]
update(R[]) <- R[i] + n_IR[i]
## Individual probabilities of transition:
p_SI[] <- 1 - exp(-beta * I[i] / N[i])
p_IR <- 1 - exp(-gamma)
## Draws from binomial distributions for numbers changing between
## compartments:
n_SI[] <- rbinom(S[i], p_SI[i])
n_IR[] <- rbinom(I[i], p_IR)
## Total population size
N[] <- S[i] + I[i] + R[i]
## Initial states:
initial(S[]) <- S_ini
initial(I[]) <- I_ini
initial(R[]) <- 0
## User defined parameters - default in parentheses:
S_ini <- user(1000)
I_ini <- user(1)
beta <- user(0.2)
gamma <- user(0.1)
## Number of replicates
nsim <- user(100)
dim(N) <- nsim
dim(S) <- nsim
dim(I) <- nsim
dim(R) <- nsim
dim(p_SI) <- nsim
dim(n_SI) <- nsim
dim(n_IR) <- nsim
## <odin_model> object generator
## Public:
## initialize: function (..., user = list(...), use_dde = FALSE, unused_user_action = NULL)
## ir: function ()
## set_user: function (..., user = list(...), unused_user_action = NULL)
## initial: function (step)
## rhs: function (step, y)
## update: function (step, y)
## contents: function ()
## transform_variables: function (y)
## engine: function ()
## run: function (step, y = NULL, ..., use_names = TRUE)
## Private:
## ptr: NULL
## use_dde: NULL
## odin: NULL
## variable_order: NULL
## output_order: NULL
## n_out: NULL
## ynames: NULL
## interpolate_t: NULL
## cfuns: list
## dll: discrete.stochastic.sir.arrays74f86426
## user: I_ini S_ini beta gamma nsim
## registration: function ()
## set_initial: function (step, y, use_dde)
## update_metadata: function ()
## Parent env: <environment: namespace:discrete.stochastic.sir.arrays74f86426>
## Locked objects: TRUE
## Locked class: FALSE
## Portable: TRUE
set.seed(1)
sir_col_transp <- paste0(sir_col, "66")
x_res <- x$run(0:100)
par(mar = c(4.1, 5.1, 0.5, 0.5), las = 1)
matplot(x_res[, 1], x_res[, -1], xlab = "Time", ylab = "Number of individuals",
type = "l", col = rep(sir_col_transp, each = 100), lty = 1)
legend("left", lwd = 1, col = sir_col, legend = c("S", "I", "R"), bty = "n")
This model is a more complex version of the previous one, which we will use to illustrate the use of all distributions mentioned in the first part: Binomial, Poisson and Multinomial.
The model is contains the following compartments:
There are no birth of natural death processes in this model. Parameters are:
The model will be written as:
$$ S_{t+1} = S_t - \beta \frac{S_t (I_{R,t} + I_{D,t})}{N_t} + \omega R_t $$
$$ E_{t+1} = E_t + \beta \frac{S_t (I_{R,t} + I_{D,t})}{N_t} - \delta E_t + \epsilon $$
IR, t + 1 = IR, t + δ(1 − μ)Et − γRIR, t + ϵ
ID, t + 1 = ID, t + δμEt − γDID, t + ϵ
Rt + 1 = Rt + γRIR, t − ωRt
Dt + 1 = Dt + γDID, t
The formulation of the model in odin
is:
## Core equations for transitions between compartments:
update(S) <- S - n_SE + n_RS
update(E) <- E + n_SE - n_EI + n_import_E
update(Ir) <- Ir + n_EIr - n_IrR
update(Id) <- Id + n_EId - n_IdD
update(R) <- R + n_IrR - n_RS
update(D) <- D + n_IdD
## Individual probabilities of transition:
p_SE <- 1 - exp(-beta * I / N)
p_EI <- 1 - exp(-delta)
p_IrR <- 1 - exp(-gamma_R) # Ir to R
p_IdD <- 1 - exp(-gamma_D) # Id to d
p_RS <- 1 - exp(-omega) # R to S
## Draws from binomial distributions for numbers changing between
## compartments:
n_SE <- rbinom(S, p_SE)
n_EI <- rbinom(E, p_EI)
n_EIrId[] <- rmultinom(n_EI, p)
p[1] <- 1 - mu
p[2] <- mu
dim(p) <- 2
dim(n_EIrId) <- 2
n_EIr <- n_EIrId[1]
n_EId <- n_EIrId[2]
n_IrR <- rbinom(Ir, p_IrR)
n_IdD <- rbinom(Id, p_IdD)
n_RS <- rbinom(R, p_RS)
n_import_E <- rpois(epsilon)
## Total population size, and number of infecteds
I <- Ir + Id
N <- S + E + I + R + D
## Initial states
initial(S) <- S_ini
initial(E) <- E_ini
initial(Id) <- 0
initial(Ir) <- 0
initial(R) <- 0
initial(D) <- 0
## User defined parameters - default in parentheses:
S_ini <- user(1000) # susceptibles
E_ini <- user(1) # infected
beta <- user(0.3) # infection rate
delta <- user(0.3) # inverse incubation
gamma_R <- user(0.08) # recovery rate
gamma_D <- user(0.12) # death rate
mu <- user(0.7) # CFR
omega <- user(0.01) # rate of waning immunity
epsilon <- user(0.1) # import case rate
## <odin_model> object generator
## Public:
## initialize: function (..., user = list(...), use_dde = FALSE, unused_user_action = NULL)
## ir: function ()
## set_user: function (..., user = list(...), unused_user_action = NULL)
## initial: function (step)
## rhs: function (step, y)
## update: function (step, y)
## contents: function ()
## transform_variables: function (y)
## engine: function ()
## run: function (step, y = NULL, ..., use_names = TRUE)
## Private:
## ptr: NULL
## use_dde: NULL
## odin: NULL
## variable_order: NULL
## output_order: NULL
## n_out: NULL
## ynames: NULL
## interpolate_t: NULL
## cfuns: list
## dll: discrete.stochastic.seirdsf90d4ad8
## user: E_ini S_ini beta delta epsilon gamma_D gamma_R mu omega
## registration: function ()
## set_initial: function (step, y, use_dde)
## update_metadata: function ()
## Parent env: <environment: namespace:discrete.stochastic.seirdsf90d4ad8>
## Locked objects: TRUE
## Locked class: FALSE
## Portable: TRUE
x <- seirds_generator$new()
seirds_col <- c("#8c8cd9", "#e67300", "#d279a6", "#ff4d4d", "#999966",
"#660000")
set.seed(1)
x_res <- x$run(0:365)
par(mar = c(4.1, 5.1, 0.5, 0.5), las = 1)
matplot(x_res[, 1], x_res[, -1], xlab = "Time", ylab = "Number of individuals",
type = "l", col = seirds_col, lty = 1)
legend("left", lwd = 1, col = seirds_col,
legend = c("S", "E", "Ir", "Id", "R", "D"), bty = "n")
Several runs can be obtained without rewriting the model, for instance, to get 100 replicates:
## [1] 366 600
## S.1 E.1 Id.1 Ir.1 R.1 D.1 S.2 E.2 Id.2 Ir.2
## 1 1000 1 0 0 0 0 1000 1 0 0
## 2 1000 1 0 0 0 0 1000 1 0 0
## 3 1000 0 0 1 0 0 1000 2 0 0
## 4 1000 0 0 1 0 0 1000 1 1 0
## 5 1000 0 0 1 0 0 1000 1 1 0
## 6 1000 0 0 1 0 0 1000 0 2 0
seirds_col_transp <- paste0(seirds_col, "1A")
par(mar = c(4.1, 5.1, 0.5, 0.5), las = 1)
matplot(0:365, x_res, xlab = "Time", ylab = "Number of individuals",
type = "l", col = rep(seirds_col_transp, 100), lty = 1)
legend("right", lwd = 1, col = seirds_col,
legend = c("S", "E", "Ir", "Id", "R", "D"), bty = "n")
It is then possible to explore the behaviour of the model using a simple function:
check_model <- function(n = 50, t = 0:365, alpha = 0.2, ...,
legend_pos = "topright") {
model <- seirds_generator$new(...)
col <- paste0(seirds_col, "33")
res <- as.data.frame(replicate(n, model$run(t)[, -1]))
opar <- par(no.readonly = TRUE)
on.exit(par(opar))
par(mar = c(4.1, 5.1, 0.5, 0.5), las = 1)
matplot(t, res, xlab = "Time", ylab = "", type = "l",
col = rep(col, n), lty = 1)
mtext("Number of individuals", side = 2, line = 3.5, las = 3, cex = 1.2)
legend(legend_pos, lwd = 1, col = seirds_col,
legend = c("S", "E", "Ir", "Id", "R", "D"), bty = "n")
}
This is a sanity check with a null infection rate and no imported case:
Another easy case: no importation, no waning immunity:
A more nuanced case: persistence of the disease with limited import, waning immunity, low severity, larger population: