R coding hitung R0 dengan metode SIR

Apa itu R0?
Grafik covid19 akan melandai atau mencapai puncaknya, jika factor R0 mendekati satu atau kurang dari satu. Sebagai contoh mudahnya jika factor R0 =3, artinya satu orang pertama akan menulari 3 orang lain, dan seterusnya berantai.  Begitu juga jika R0=2, nah jika R0 lebih kecil dari satu atau negative  itu artinya covid19 berhenti.

Data
Diambil data covid19 Indonesia, antara  "2020-03-02"-"2020-05-02"

Coding R
Sebelum menjalankan Coding R, install dulu Rstudio

library yang diperlukan

library(coronavirus)
library(tidyr)
library(tableHTML)
library(coronavirus)
library(dplyr)
library(magrittr)
library(lubridate)
library(deSolve)

Ambil data covid19

data(coronavirus)
tt <- head(coronavirus)
head(tt)

##         date province     country lat long      type cases
## 1 2020-01-22          Afghanistan  33   65 confirmed     0
## 2 2020-01-23          Afghanistan  33   65 confirmed     0
## 3 2020-01-24          Afghanistan  33   65 confirmed     0
## 4 2020-01-25          Afghanistan  33   65 confirmed     0
## 5 2020-01-26          Afghanistan  33   65 confirmed     0
## 6 2020-01-27          Afghanistan  33   65 confirmed     0

Rumus dasar SIR

SIR <- function(time, state, parameters) {
  par <- as.list(c(state, parameters))
  with(par, {
    dS <- -beta * I * S / N
    dI <- beta * I * S / N - gamma * I
    dR <- gamma * I
    list(c(dS, dI, dR))
  })
}

Proses perhitungan covid19 dg metode SIR

#library(dplyr)
#library(magrittr)

df <- coronavirus %>%
  filter(country == "Indonesia") %>%
  group_by(date, type) %>%
  summarise(total = sum(cases, na.rm = TRUE)) %>%
  pivot_wider(
    names_from = type,
    values_from = total
  ) %>%
  arrange(date) %>%
  ungroup() %>%
  mutate(active = confirmed - death - recovered) %>%
  mutate(
    confirmed_cum = cumsum(confirmed),
    death_cum = cumsum(death),
    recovered_cum = cumsum(recovered),
    active_cum = cumsum(active)
  )
td <- tail(df,50)
td

## # A tibble: 50 x 9
##    date       confirmed death recovered active confirmed_cum death_cum
##    <date>         <int> <int>     <int>  <int>         <int>     <int>
##  1 2020-04-05       181     7        14    160          2273       198
##  2 2020-04-06       218    11        28    179          2491       209
##  3 2020-04-07       247    12        12    223          2738       221
##  4 2020-04-08       218    19        18    181          2956       240
##  5 2020-04-09       337    40        30    267          3293       280
##  6 2020-04-10       219    26        30    163          3512       306
##  7 2020-04-11       330    21         4    305          3842       327
##  8 2020-04-12       399    46        73    280          4241       373
##  9 2020-04-13       316    26        21    269          4557       399
## 10 2020-04-14       282    60        46    176          4839       459
## # ... with 40 more rows, and 2 more variables: recovered_cum <int>,
## #   active_cum <int>

Menentukan start data dan end data

#library(lubridate)

sir_start_date <- "2020-03-02"
sir_end_date <- "2020-05-02"

Infected <- subset(df, date >= ymd(sir_start_date) & date <= ymd(sir_end_date))$active_cum
# Create an incrementing Day vector the same length as our
# cases vector
Day <- 1:(length(Infected))

Menentukan fungsi awal S, I, R masukkan nilai N=jumlah penduduk Indonesia

N <- 268074600 #Jumlah penduduk Indonesia, https://id.wikipedia.org/wiki/
init <- c(
  S = N - Infected[1],
  I = Infected[1],
  R = 0
)

RSS <- function(parameters) {
  names(parameters) <- c("beta", "gamma")
  out <- ode(y = init, times = Day, func = SIR, parms = parameters)
  fit <- out[, 3]
  sum((Infected - fit)^2)
  }

#library(deSolve)

Opt <- optim(c(0.5, 0.5),
  RSS,
  method = "L-BFGS-B",
  lower = c(0, 0),
  upper = c(1, 1)
)

Check for convergence

Opt$message

## [1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"

#[1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"
#---Convergence is confirmed. Now we can examine the fitted values for
Opt_par <- setNames(Opt$par, c("beta", "gamma"))
Opt_par

##      beta     gamma 
## 0.5715389 0.4284613

Hasil perhitungan RO berdasarkan data covid19 Indonesia, antara 2020-03-02 sampai 2020-05-02

R0 <- as.numeric(Opt_par[1] / Opt_par[2])
R0

## [1] 1.333934