2  SEIR

2.1 Code

library(odin)
library(tidyr)
library(ggplot2)
library(RColorBrewer)
library(patchwork)

Let reuse the deterministic SIR model code from Listing 1.1 and write code for the deterministic SEIR model.

Listing 2.1: Deterministic SEIR model

odin_seir <- odin::odin({
  # Derivatives
  deriv(S) <- -beta * S * I
  deriv(E) <- beta * S * I - sigma * E
  deriv(I) <- sigma * E - gamma * I
  deriv(R) <- gamma * I
  deriv(CInc) <- sigma * E
  
  # Initial conditions
  initial(S) <- S_init
  initial(E) <- E_init
  initial(I) <- I_init
  initial(R) <- R_init
  initial(CInc) <- CInc_init
  
  # Parameters and initial values
  beta <- user(0.004)
  sigma <- user(0.5)
  gamma <- user(0.5)
  S_init <- user(369)
  E_init <- user(0)
  I_init <- user(1)
  R_init <- user(0)
  CInc_init <- user(0)
})

seir_mod <- function(beta, sigma, gamma, S0, E0, I0, R0, times) {
  # Set values for a new run
  odin_run <- odin_seir$new(beta = beta, sigma = sigma, gamma = gamma, S_init = S0, E_init = E0, I_init = I0, R_init = R0)
  
  # Run the model
  out <- data.frame(odin_run$run(times))

  # Compute incidence from cumulative incidence
  out$Inc <- c(I0, diff(out$CInc))
  out$CInc <- NULL
  out
}

pred1 <- sir_mod(beta = 0.004, gamma = 0.5, S0 = 999, I0 = 1, R0 = 0, times = 1:20)
pred2 <- seir_mod(beta = 0.004, sigma = 0.5, gamma = 0.5, S0 = 999, E0 = 0, I0 = 1, R0 = 0, times = 1:20)

df_plot <- pivot_longer(pred1, cols = S:Inc, names_to = "comp", values_to = "n")

my_palette <- brewer.pal(11, "PuOr")[c(10, 4, 3, 8)]

p1 <- ggplot(df_plot, aes(x = t, y = n, color = comp)) +
  geom_line(linewidth = 1.2) +
  scale_color_manual(values = my_palette, breaks = c("S", "Inc", "I", "R")) +
  labs(color = NULL, y = NULL, x = NULL) +
  theme_minimal() +
  theme(legend.position = "none")

df_plot <- pivot_longer(pred2, cols = S:Inc, names_to = "comp", values_to = "n")

my_palette <- brewer.pal(11, "PuOr")[c(10, 1, 4, 3, 8)]

p2 <- ggplot(df_plot, aes(x = t, y = n, color = comp)) +
  geom_line(linewidth = 1.2) +
  scale_color_manual(values = my_palette, breaks = c("S", "E", "Inc", "I", "R")) +
  labs(color = NULL, y = NULL, x = "Time") +
  theme_minimal() +
  theme(legend.position = "bottom")

p1 / p2