Nothing
odin discrete models
In odin: ODE Generation and Integration
{r setup, echo = FALSE, results = "hide"}
lang_output <- function(x, lang) {
cat(c(sprintf("%s", lang), x, "```"), sep = "\n")
}
c_output <- function(x) lang_output(x, "cc")
r_output <- function(x) lang_output(x, "r")
plain_output <- function(x) lang_output(x, "plain")
knitr::opts_chunk$set(
fig.width = 7,
fig.height = 5)
options(odin.verbose = FALSE)
# Discrete compartmental models in a nutshell
## From continuous to discrete time
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 $\beta$ is an infection rate and $\gamma$ 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
$$
$$
R_{t+1} = R_t + \gamma I_t
$$
## Stochastic processes
The discrete model above remains deterministic: for given values of
the rates $\beta$ and $\gamma$, 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.
### Binomial distribution
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:
1. 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 $\lambda$, then
the corresponding probability of this individual leaving the
compartment in one time unit is $p = 1 - e^{- \lambda}$.
2. 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 \rightarrow 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 distribution
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
$\epsilon$. 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 $\mathcal{P}(\epsilon)$
- individual move from $I$ to $R$ according to a Binomial
distribution $\mathcal{B}(I_t, 1 - e^{- \gamma})$
### Multinomial distribution
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 $\lambda_A$, $\lambda_B$ and $\lambda_C$. The
workflow to handle these transitions will be:
1. determine the total number of individuals leaving $X$; this is
done by summing the rates ($\lambda = \lambda_A + \lambda_B +
\lambda_C$) to compute the *per capita* probability of leaving $X$
$(p_(X \rightarrow ...) = 1 - e^{- \lambda})$, and drawing the
number of individuals leaving $X$ ($n_{_(X \rightarrow ...)}$) from
a binomial distribution $n_{(X \rightarrow ...)} \sim B(X, p_(X
\rightarrow ...))$
2. 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}$
3. determine the numbers of individuals moving to $A$, $B$ and $C$
using a Multinomial distribution: $n_{(X \rightarrow A, B, C)} \sim
\mathcal{M}(n_{(X \rightarrow ...)}, p_A, p_B, p_C)$
# Implementation using `odin`
## Deterministic SIR model
We start by loading the `odin` code for a discrete, stochastic SIR
model:
``` {r load_sir}
path_sir_model <- system.file("examples/discrete_deterministic_sir.R", package = "odin")
``` {r echo = FALSE, results = "asis"}
r_output(readLines(path_sir_model))
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:
``` {r }
sir_generator <- odin::odin(path_sir_model)
sir_generator
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_generatoris an R6 generator that can be used to create an instance of the model:
generate an instance of the model:
x <- sir_generator$new()
x
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.:
``` {r sir-deterministic, fig.cap = "An example of deterministic, discrete-time SIR model
"}
sir_col <- c("#8c8cd9", "#cc0044", "#999966")
x$run(0:10)
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")
## Stochastic SIR model
The stochastic equivalent of the previous model can be formulated
in `odin` as follows:
``` {r load_sir_s}
path_sir_model_s <- system.file("examples/discrete_stochastic_sir.R", package = "odin")
``` {r echo = FALSE, results = "asis"}
r_output(readLines(path_sir_model_s))
We can use the same workflow as before to run the model, using 10
initial infected individuals (`I_ini = 10`):
``` {r }
sir_s_generator <- odin::odin(path_sir_model_s)
sir_s_generator
x <- sir_s_generator$new(I_ini = 10)
``` {r sir-stochastic_1, fig.cap = "An example of stochastic, discrete-time SIR model
"}
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:
``` {r }
path_sir_model_s_a <- system.file("examples/discrete_stochastic_sir_arrays.R", package = "odin")
``` {r echo = FALSE, results = "asis"}
r_output(readLines(path_sir_model_s_a))
``` {r echo = TRUE}
sir_s_a_generator <- odin::odin(path_sir_model_s_a)
sir_s_a_generator
x <- sir_s_a_generator$new()
``` {r sir-stochastic_100, fig.cap = "100 replicates of a stochastic, discrete-time SIR model
"}
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")
## A stochastic SEIRDS model
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:
- $S$: susceptibles
- $E$: exposed, i.e. infected but not yet contagious
- $I_R$: infectious who will survive
- $I_D$: infectious who will die
- $R$: recovered
- $D$: dead
There are no birth of natural death processes in this model. Parameters are:
- $\beta$: rate of infection
- $\delta$: rate at which symptoms appear (i.e inverse of mean incubation
period)
- $\gamma_R$: recovery rate
- $\gamma_D$: death rate
- $\mu$: case fatality ratio (proportion of cases who die)
- $\epsilon$: import rate of infected individuals (applies to $E$ and $I$)
- $\omega$: rate waning immunity
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
$$
$$
I_{R,t+1} = I_{R,t} + \delta (1 - \mu) E_t - \gamma_R I_{R,t} + \epsilon
$$
$$
I_{D,t+1} = I_{D,t} + \delta \mu E_t - \gamma_D I_{D,t} + \epsilon
$$
$$
R_{t+1} = R_t + \gamma_R I_{R,t} - \omega R_t
$$
$$
D_{t+1} = D_t + \gamma_D I_{D,t}
$$
The formulation of the model in `odin` is:
``` {r load_seirds}
path_seirds_model <- system.file("examples/discrete_stochastic_seirds.R", package = "odin")
``` {r echo = FALSE, results = "asis"}
r_output(readLines(path_seirds_model))
``` {r seirds}
seirds_generator <- odin::odin(path_seirds_model)
seirds_generator
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:
``` {r seirds_100, fig.cap = "100 replicates of a stochastic, discrete-time SEIRDS model
"}
x_res <- as.data.frame(replicate(100, x$run(0:365)[, -1]))
dim(x_res)
x_res[1:6, 1:10]
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:
``` {r custom_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:
``` {r fig.cap = "Stochastic SEIRDS model: sanity check with no infections
"}
check_model(beta = 0, epsilon = 0)
Another easy case: no importation, no waning immunity:
``` {r fig.cap = "<i>Stochastic SEIRDS model: no importation or waning immunity</i><br>"}
check_model(epsilon = 0, omega = 0)
A more nuanced case: persistence of the disease with limited
import, waning immunity, low severity, larger population:
{r fig.cap = "<i>Stochastic SEIRDS model: endemic state in a larger population</i><br>"}
check_model(t = 0:(365 * 3), epsilon = 0.1, beta = .2, omega = .01,
mu = 0.005, S_ini = 1e5)
Try the odin package in your browser
Any scripts or data that you put into this service are public.
odin documentation built on Oct. 2, 2023, 5:07 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
{r setup, echo = FALSE, results = "hide"}
lang_output <- function(x, lang) {
cat(c(sprintf("%s", lang), x, "```"), sep = "\n")
}
c_output <- function(x) lang_output(x, "cc")
r_output <- function(x) lang_output(x, "r")
plain_output <- function(x) lang_output(x, "plain")
knitr::opts_chunk$set(
fig.width = 7,
fig.height = 5)
options(odin.verbose = FALSE)
# Discrete compartmental models in a nutshell ## From continuous to discrete time 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 $\beta$ is an infection rate and $\gamma$ 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 $$ $$ R_{t+1} = R_t + \gamma I_t $$ ## Stochastic processes The discrete model above remains deterministic: for given values of the rates $\beta$ and $\gamma$, 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. ### Binomial distribution 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: 1. 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 $\lambda$, then the corresponding probability of this individual leaving the compartment in one time unit is $p = 1 - e^{- \lambda}$. 2. 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 \rightarrow 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 distribution 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 $\epsilon$. 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 $\mathcal{P}(\epsilon)$ - individual move from $I$ to $R$ according to a Binomial distribution $\mathcal{B}(I_t, 1 - e^{- \gamma})$ ### Multinomial distribution 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 $\lambda_A$, $\lambda_B$ and $\lambda_C$. The workflow to handle these transitions will be: 1. determine the total number of individuals leaving $X$; this is done by summing the rates ($\lambda = \lambda_A + \lambda_B + \lambda_C$) to compute the *per capita* probability of leaving $X$ $(p_(X \rightarrow ...) = 1 - e^{- \lambda})$, and drawing the number of individuals leaving $X$ ($n_{_(X \rightarrow ...)}$) from a binomial distribution $n_{(X \rightarrow ...)} \sim B(X, p_(X \rightarrow ...))$ 2. 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}$ 3. determine the numbers of individuals moving to $A$, $B$ and $C$ using a Multinomial distribution: $n_{(X \rightarrow A, B, C)} \sim \mathcal{M}(n_{(X \rightarrow ...)}, p_A, p_B, p_C)$ # Implementation using `odin` ## Deterministic SIR model We start by loading the `odin` code for a discrete, stochastic SIR model: ``` {r load_sir} path_sir_model <- system.file("examples/discrete_deterministic_sir.R", package = "odin")
``` {r echo = FALSE, results = "asis"} r_output(readLines(path_sir_model))
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: ``` {r } sir_generator <- odin::odin(path_sir_model) sir_generator
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_generatoris an R6 generator that can be used to create an instance of the model:
generate an instance of the model:
x <- sir_generator$new() x
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.:
``` {r sir-deterministic, fig.cap = "An example of deterministic, discrete-time SIR model
"}
sir_col <- c("#8c8cd9", "#cc0044", "#999966")
x$run(0:10)
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")
## Stochastic SIR model The stochastic equivalent of the previous model can be formulated in `odin` as follows: ``` {r load_sir_s} path_sir_model_s <- system.file("examples/discrete_stochastic_sir.R", package = "odin")
``` {r echo = FALSE, results = "asis"} r_output(readLines(path_sir_model_s))
We can use the same workflow as before to run the model, using 10 initial infected individuals (`I_ini = 10`): ``` {r } sir_s_generator <- odin::odin(path_sir_model_s) sir_s_generator x <- sir_s_generator$new(I_ini = 10)
``` {r sir-stochastic_1, fig.cap = "An example of stochastic, discrete-time SIR model
"}
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: ``` {r } path_sir_model_s_a <- system.file("examples/discrete_stochastic_sir_arrays.R", package = "odin")
``` {r echo = FALSE, results = "asis"} r_output(readLines(path_sir_model_s_a))
``` {r echo = TRUE}
sir_s_a_generator <- odin::odin(path_sir_model_s_a)
sir_s_a_generator
x <- sir_s_a_generator$new()
``` {r sir-stochastic_100, fig.cap = "100 replicates of a stochastic, discrete-time SIR model
"}
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")
## A stochastic SEIRDS model 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: - $S$: susceptibles - $E$: exposed, i.e. infected but not yet contagious - $I_R$: infectious who will survive - $I_D$: infectious who will die - $R$: recovered - $D$: dead There are no birth of natural death processes in this model. Parameters are: - $\beta$: rate of infection - $\delta$: rate at which symptoms appear (i.e inverse of mean incubation period) - $\gamma_R$: recovery rate - $\gamma_D$: death rate - $\mu$: case fatality ratio (proportion of cases who die) - $\epsilon$: import rate of infected individuals (applies to $E$ and $I$) - $\omega$: rate waning immunity 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 $$ $$ I_{R,t+1} = I_{R,t} + \delta (1 - \mu) E_t - \gamma_R I_{R,t} + \epsilon $$ $$ I_{D,t+1} = I_{D,t} + \delta \mu E_t - \gamma_D I_{D,t} + \epsilon $$ $$ R_{t+1} = R_t + \gamma_R I_{R,t} - \omega R_t $$ $$ D_{t+1} = D_t + \gamma_D I_{D,t} $$ The formulation of the model in `odin` is: ``` {r load_seirds} path_seirds_model <- system.file("examples/discrete_stochastic_seirds.R", package = "odin")
``` {r echo = FALSE, results = "asis"} r_output(readLines(path_seirds_model))
``` {r seirds}
seirds_generator <- odin::odin(path_seirds_model)
seirds_generator
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:
``` {r seirds_100, fig.cap = "100 replicates of a stochastic, discrete-time SEIRDS model
"}
x_res <- as.data.frame(replicate(100, x$run(0:365)[, -1]))
dim(x_res)
x_res[1:6, 1:10]
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: ``` {r custom_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:
``` {r fig.cap = "Stochastic SEIRDS model: sanity check with no infections
"}
check_model(beta = 0, epsilon = 0)
Another easy case: no importation, no waning immunity:
``` {r fig.cap = "<i>Stochastic SEIRDS model: no importation or waning immunity</i><br>"}
check_model(epsilon = 0, omega = 0)
A more nuanced case: persistence of the disease with limited import, waning immunity, low severity, larger population:
{r fig.cap = "<i>Stochastic SEIRDS model: endemic state in a larger population</i><br>"}
check_model(t = 0:(365 * 3), epsilon = 0.1, beta = .2, omega = .01,
mu = 0.005, S_ini = 1e5)
Try the odin package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.