Kapan Puncak Covid19 di Indonesia ?

Prediksi Covid19, kapan mencapai puncaknya?

Selamat Idhul Fitri from Home, menjawab pertanyaan kapan covid19 mencapai puncaknya? Atau kapan grafiknya melandai? Kalau ada pejabat yang mengatakan, kita sudah mulai lega karena R0(RNol) sudah mendekati 1(satu). Apa artinya? Dan benarkah itu?

Mbak Mona waspada Covid19

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.

Grafik Prediksi Covid19, 150 hr ke depan

R0 dihitung dari turunan Sucestible(ODP), Infected(Terinfeksicovid19), Recovery(Sembuh)  perhitungan ini dalam data science dikenal metode SIR dengan R coding. Detail coding ada di bawah

Apa benar saat ini di Indonesia R0(RNol) = 1?

Saya menghitung dengan metode SIR, hasilnya 1.3339341 selanjutnya saya membuat prediksi 150 hari kedepan  dan hasinya Covid19 akan mencapai puncaknya pada akhir bulan Juli 2020. Analisa ini berdasar real data Covid19 diambil dari Library(coronavirus) data ini sama dengan data yang ditayangkan situs covid19.co.id.Dengan acuan jumlah penduduk Indonesia saat ini 267juta(wikipedia.com). Berikut saya tuliskan script R nya:

ui

library(shiny)
library(miniUI)
library(coronavirus)
library(magrittr)
library(tidyr)
library(dplyr)
library(tableHTML)
library(shinyWidgets)
  ui <-
   miniPage(
 #           setBackgroundColor(
  #  color = c("#F7FBFF", "#2171B5"),
  #  gradient = "radial",
  #  direction = c("top", "left")
  #),
   # gadgetTitleBar("VHD Analytic"),
    gadgetTitleBar("___PREDIKSI COVID-19 INDONESIA 150 HARI KE DEPAN___",setBackgroundColor(
    color = c("blue", "magenta"),
    gradient = "radial",
    direction = c("top", "left"))),
 
    miniTabstripPanel(
    miniTabPanel("Data Global Covid19", icon = icon("table"),
      miniContentPanel(       
       verbatimTextOutput("tt"))),
 
    miniTabPanel("Data Covid19 Indonesia", icon = icon("table"),
      miniContentPanel(       
       verbatimTextOutput("td"))),
     miniTabPanel("Metode SIR", icon = icon("list-alt"),
      miniContentPanel(       
       dataTableOutput("About"),
               h3("Menghitung R0"),
 p("Illustrasi Metode SIR",align = 'Justify'),   
 p("SIR <- function(time, state, parameters) {",align = 'Justify'),
 p("par <- as.list(c(state, parameters))",align = 'Justify'),
 p("with(par, {",align = 'Justify'),
 p("dS <- -beta * I * S / N",align = 'Justify'),
 p("dI <- beta * I * S / N - gamma * I",align = 'Justify'),
 p("dR <- gamma * I",align = 'Justify'),
 p("list(c(dS, dI, dR))",align = 'Justify'),
 p("})",align = 'Justify'),
 p("...dst didapat beta 0.5715389 gamma 0.4284613",align = 'Justify'),
 p("...dst beta/gamma didapat R0 0.333934",align = 'Justify'),
 h3("Prediksi 150 hr kedepan",align = 'Justify'),
 p("Data Covid19 Indonesia, sumber library(coronavirus), datanya sama dg covid19.co.id",align = 'Justify'),
 p("Jumlah penduduk Indonesia 268074600, sumber: wikipedia.org",align = 'Justify'))),
 miniTabPanel("Variable R0",icon = icon("area-chart"),
       miniContentPanel(       
        verbatimTextOutput(("R0")))),
     miniTabPanel("Grafik covid19",icon = icon("area-chart"),
    miniContentPanel(
       plotOutput("p1"))),
miniTabPanel("Prediksi 150 hr",icon = icon("area-chart"),
    miniContentPanel(
       plotOutput("p2")))
  ))

global

library(shiny)
library(miniUI)
library(coronavirus)
library(magrittr)
library(tidyr)
library(tableHTML)
library(dplyr)
library(shinyWidgets)
library(coronavirus)
data(coronavirus)
tt <- head(coronavirus,50)
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))
  })
}
#write.csv(coronavirus,"d:/data/coronavirus.csv")
library(dplyr)
library(magrittr)
# extract the cumulative incidence
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)
## A tibble: 10 x 9
##   date       confirmed death recovered active confirmed_cum death_cum recovered_cum ##active_cum
##   <date>         <int> <int>     <int>  <int>         <int>     <int>         <int>      <int>
## 1 2020-05-15       490    33       285    172         16496      1076          3803      11617
## 2 2020-05-16       529    13       108    408         17025      1089          3911      12025
## 3 2020-05-17       489    59       218    212         17514      1148          4129      12237
## 4 2020-05-18       496    43       195    258         18010      1191          4324      12495
## 5 2020-05-19       486    30       143    313         18496      1221          4467      12808

# put the daily cumulative incidence numbers for Belgium from
# Feb 4 to March 30 into a vector called Infected
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))

# now specify initial values for N, S, I and R
N <- 268074600 #https://id.wikipedia.org/wiki/Daftar_negara_menurut_jumlah_penduduk
init <- c(
  S = N - Infected[1],
  I = Infected[1],
  R = 0
)
# optimised for the best fit to the incidence data
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)
}

# install.packages("deSolve")
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"
#---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

R0 <- as.numeric(Opt_par[1] / Opt_par[2])
R0
##[1] 1.333934 faktor R0 hasil perhitungan dengan metode SIR

# time in days for predictions
t <- 1:as.integer(ymd(sir_end_date) + 1 - ymd(sir_start_date))

# get the fitted values from our SIR model
fitted_cumulative_incidence <- data.frame(ode(
  y = init, times = t,
  func = SIR, parms = Opt_par
))

# add a Date column and the observed incidence data
library(dplyr)
fitted_cumulative_incidence <- fitted_cumulative_incidence %>%
  mutate(
    Date = ymd(sir_start_date) + days(t - 1),
    Country = "Indonesia",
    cumulative_incident_cases = Infected
  )

# plot the data
library(ggplot2)
p1 <- fitted_cumulative_incidence %>%
  ggplot(aes(x = Date)) +
  geom_line(aes(y = I), colour = "red") +
  geom_point(aes(y = cumulative_incident_cases), colour = "blue") +
  labs(
    y = "Cumulative incidence",
    title = "COVID-19 Cumulative Indonesia, Indonesia",
    subtitle = "(Red = metode SIR, blue = data observasi)"
  ) +
  theme_minimal() #ok chart 1
# plot the data
t <- 1:150
fitted_cumulative_incidence <- data.frame(ode(
  y = init, times = t,
  func = SIR, parms = Opt_par
))
# add a Date column and join the observed incidence data
fitted_cumulative_incidence <- fitted_cumulative_incidence %>%
  mutate(
    Date = ymd(sir_start_date) + days(t - 1),
    Country = "Indonesia",
    cumulative_incident_cases = c(Infected, rep(NA, length(t) - length(Infected)))
  )
p2 <- fitted_cumulative_incidence %>%
  ggplot(aes(x = Date)) +
  geom_line(aes(y = I, colour = "red")) +
  geom_line(aes(y = S, colour = "black")) +
  geom_line(aes(y = R, colour = "green")) +
  geom_point(aes(y = cumulative_incident_cases, colour = "blue")) +
  scale_y_log10(labels = scales::comma) +
  labs(
    y = "Persons",
    title = "Prediksi COVID19 di Indonesia, 150 hr kedepan, puncaknya akhir Juli 2020, ditandai garis(infectious)"
  ) +
  scale_colour_manual(
    name = "",
    values = c(red = "red", black = "black", green = "green", blue = "blue"),
    labels = c("Susceptible", "Observed", "Recovered", "Infectious")
  ) +
  theme_minimal()#ok chart 

server

server <- function(input, output, session) {
   output$tt <- renderPrint({
    tt
  }) 
  output$td <- renderPrint({
    td
  }) 
  output$R0<- renderPrint({
    R0
  })
    output$p1 <- renderPlot({
    plot(p1)
  })
   output$p2 <- renderPlot({
    plot(p2)
  })
    observeEvent(input$done, {
    stopApp(TRUE)
  })
}

Klik:
https://bambangpe.shinyapps.io/covid19-prediction/
https://bambangpe.shinyapps.io/kursus-r2/
https://covid19.go.id/