source:https://epiverse-trace.github.io/tutorials-late/index.html

I. Summary and Setup

library(tidyverse)   # Sudah mencakup dplyr, readr, ggplot2, dll.
library(socialmixr)
library(epidemics)
library(scales)

II. Load Dataset

ex1 <- read_csv("c:/users/sonk/Documents/Kornelius Son/Sintax R/sarscov2.csv")
ex2 <- read_csv("c:/users/sonk/Documents/Kornelius Son/Sintax R/ebola_cases.csv")

III. Simulating Disease Spread

Let us assume that people are, on average, infectious for 8 days. Using this assumption, we can calculate the recovery rate (γ), as follows: with γ is equal to 1/Duration of the infectiousness

gamma <- 1 / 8  # Recovery rate, where infectious period is 8 days
gamma
## [1] 0.125

IV. Load contact and population data

Using the socialmixr package, we create the contact matrix for the UK for the age bins: 0-20 yo, 20-40 yo, and >=40 yo.

# Load the polymod data and create a contact matrix
polymod <- socialmixr::polymod

# Generate the contact matrix for UK with specified age limits
contact_data <- socialmixr::contact_matrix(
     survey = polymod,
     countries = "United Kingdom",
     age.limits = c(0, 20, 40),
     symmetric = TRUE
)
#preparing contact matrix
contact_matrix <-t(contact_data$matrix)
contact_matrix
##                  age.group
## contact.age.group   [0,20)  [20,40)      40+
##           [0,20)  7.883663 2.794154 1.565665
##           [20,40) 3.120220 4.854839 2.624868
##           40+     3.063895 4.599893 5.005571

V. Initial Condition

the initial condition are the proportion of individuals in each state S,E,I,R. In this example we have three groups age, let assume that in the youngest age category, 0ne millio individuals are infectious, and the remaining age groups are infection free. S(0)= 1-1/1,000,000 ; E(0)=0, I(0)=0, R(0)=0

initial_i <- 1e-6
initial_conditions_inf <- c(
     S = 1 - initial_i, E = 0, I = initial_i, R = 0, V = 0
)
initial_conditions_free <- c(
     S = 1, E = 0, I = 0, R = 0, V = 0
)                                
#combine the initial conditions vector into one matrix
initial_conditions <-rbind(
     initial_conditions_inf, #age group 1
     initial_conditions_free, #age group 2
     initial_conditions_free #age group 3
)

rownames(initial_conditions) <-rownames(contact_matrix)
initial_conditions
##                S E     I R V
## [0,20)  0.999999 0 1e-06 0 0
## [20,40) 1.000000 0 0e+00 0 0
## 40+     1.000000 0 0e+00 0 0
#Population structure
demography_vector <- contact_data$demography$population
names(demography_vector) <- rownames(contact_matrix)
demography_vector
##   [0,20)  [20,40)      40+ 
## 14799290 16526302 28961159
uk_population <-population(
     name = "UK",
     contact_matrix = contact_matrix,
     demography_vector = demography_vector,
     initial_conditions = initial_conditions
)

VI. Model parameters

To run our model we need to specify the parameter:

We will simulate a strain of influenza with pandemic potential with \(R_0\)$=1.46, with a pre-infectious period of 3 days and infectious period of 7 days. Therefore our inputs will be:

# time periods
preinfectious_period <- 3.0
infectious_period <- 7.0
basic_reproduction <- 1.46

# rates
infectiousness_rate <- 1.0 / preinfectious_period
recovery_rate <- 1.0 / infectious_period
transmission_rate <- basic_reproduction / infectious_period

Now we are ready to run our model using model_default() from the epidemics package. Let’s specify time_end=600 to run the model for 600 days.

VII. Running Model

output <- model_default(
     # population
     population = uk_population,
     # rates
     transmission_rate = transmission_rate,
     infectiousness_rate = infectiousness_rate,
     recovery_rate = recovery_rate,
     # time
     time_end = 600, increment = 1.0
)
head(output)
##     time demography_group compartment    value
##    <num>           <char>      <char>    <num>
## 1:     0           [0,20) susceptible 14799275
## 2:     0          [20,40) susceptible 16526302
## 3:     0              40+ susceptible 28961159
## 4:     0           [0,20)     exposed        0
## 5:     0          [20,40)     exposed        0
## 6:     0              40+     exposed        0

Note: This model also has the capability to incorporate vaccination and tracks the number of vaccinated individuals over time. Although we haven’t included any vaccination in this setup, there is still a vaccinated compartment in the output, which remains empty. The use of vaccination will be discussed in future tutorials.

The output of our model shows the number of individuals in each compartment for each age group over time. We can visualize the number of infectious individuals (those in the I class) over time.

output %>%
     filter(compartment == "infectious") %>%
     ggplot() +
     geom_line(
          aes(
               x = time,
               y = value,
               colour = demography_group
          )
     ) +
     scale_y_continuous(
          labels = scales::comma
     ) +
     theme_bw() +
     labs(
          x = "Simulation time (days)",
          linetype = "Compartment",
          y = "Individuals"
     )