Model SIR adalah model epidemik yang paling umum digunakan dalam melihat persebaran penyakit. Model SIR yang diperkenalkan oleh Kermack dan McKendrick pada tahun 1927 dan pertama kali digunakan untuk melihat dinamika penyebaran populasi akibat penyakit menular.

Membaca data KasusHarian.csv

data.corona <- read.csv("KasusHarian.csv")
head(data.corona)
##                 ï..DT Kasus..Kumulatif. Kasus.Baru Sembuh..baru.
## 1 2020-03-02 00:00:00                 2          2            NA
## 2 2020-03-03 00:00:00                 2         NA            NA
## 3 2020-03-04 00:00:00                 2         NA            NA
## 4 2020-03-05 00:00:00                 2         NA            NA
## 5 2020-03-06 00:00:00                 4          2            NA
## 6 2020-03-07 00:00:00                 4         NA            NA
##   Meninggal..baru.
## 1               NA
## 2               NA
## 3               NA
## 4               NA
## 5               NA
## 6               NA

Membuat data harian dan menggabungkannya dalam data.corona

hari <- 1:26
data.corona <- cbind(data.corona,hari)

Mengubah data NA menjadi bernilai 0

data.corona[is.na(data.corona)] <- 0

Membuat data kumulatif pada kolom sembuh.baru

data.corona$sembuh.kumulatif <- cumsum(data.corona$Sembuh..baru.)

Visualisasi kasus.kumulatif

library(ggplot2)
## Warning: package 'ggplot2' was built under R version 3.5.3
plot.corona <- ggplot(data = data.corona, aes(x=hari,y=Kasus..Kumulatif.))+
  geom_point()+geom_smooth(method = "loess", formula = y~x,se=T)+ylab("Jumlah")+
  labs(title = "Total Infeksi COVID-19 Indonesia")+theme(axis.text.x=element_text(
    angle = 45,hjust = 1
  ))
plot.corona

Membuat pemodelan penyebaran virus corona dengan model SIR berdasarkan jumlah penduduk Indonesia Mengubah data ke tipe numerik.

Infected <- as.numeric(data.corona$Kasus..Kumulatif.)
Day <- as.numeric(data.corona$hari)
N <- 238518000
R <-as.numeric(data.corona$sembuh.kumulatif)

Model SIR

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

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
}

opt <- optim(c(0.5,0.5),RSS,method="L-BFGS-B",lower=c(0,0),upper = c(1,1))
opt$message
## [1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"

Memanggil nilai beta dan gamma

opt_par <- setNames(opt$par,c("beta","gamma"))
opt_par
##      beta     gamma 
## 0.6329087 0.3670913

Membuat data persebaran dengan waktu persebaran 90 hari

t <- 1:90
fit <- data.frame(ode(y=init, times=t,func=SIR, parms = opt_par))

Visualisasi model SIR

col <- 1:3
matplot(fit$time, fit[,2:4],type="l",xlab = "Day", ylab="Jumlah Manusia",lwd = 2,lty=1,
        col=col,log="y")
## Warning in xy.coords(x, y, xlabel, ylabel, log = log): 1 y value <= 0
## omitted from logarithmic plot
points(Day,Infected)
legend("bottomright",c("Rentan","Terinfeksi","Sembuh"),lty=1,lwd=2,col=col)
title("Model SIR COVID-19 Indonesia",outer=T,line=-2)

Memanggil nilai RO

RO <- setNames(opt_par["beta"]/opt_par["gamma"],"RO")
print(RO)
##       RO 
## 1.724118

Menetukan puncak persebaran dan waktunya

fit[fit$I==max(fit$I),"I",drop=F]
##           I
## 69 24788899

Menentukan jumlah kematian degan tingkat rasio kematian 0.4%

max_infected <- max(fit$I)
max_infected*0.004 
## [1] 99155.6

Membuat pemodelan virus corona dengan model SIR berdasarkan jumlah penduduk di pulau jawa. Mengubah data ke tipe numerik

Infected <- as.numeric(data.corona$Kasus..Kumulatif.)
Day <- as.numeric(data.corona$hari)
N <- 140000000
R <-as.numeric(data.corona$sembuh.kumulatif)

Model SIR

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

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
}

opt <- optim(c(0.5,0.5),RSS,method="L-BFGS-B",lower=c(0,0),upper = c(1,1))
opt$message
## [1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"

Memanggil nilai beta dan gamma

opt_par <- setNames(opt$par,c("beta","gamma"))
opt_par
##     beta    gamma 
## 0.632909 0.367091

Membuat data persebaran dengan waktu persebaran 90 hari

t <- 1:90
fit <- data.frame(ode(y=init, times=t,func=SIR, parms = opt_par))

Visualisasi model SIR

col <- 1:3
matplot(fit$time, fit[,2:4],type="l",xlab = "Day", ylab="Jumlah Manusia",lwd = 2,lty=1,
        col=col,log="y")
## Warning in xy.coords(x, y, xlabel, ylabel, log = log): 1 y value <= 0
## omitted from logarithmic plot
points(Day,Infected)
legend("bottomright",c("Rentan","Terinfeksi","Sembuh"),lty=1,lwd=2,col=col)
title("Model SIR COVID-19 Pulau Jawa",outer=T,line=-2)

Memanggil nilai RO

RO <- setNames(opt_par["beta"]/opt_par["gamma"],"RO")
print(RO)
##      RO 
## 1.72412

Menetukan puncak persebaran dan waktunya

fit[fit$I==max(fit$I),"I",drop=F]
##           I
## 67 14549598

Menentukan jumlah kematian degan tingkat rasio kematian 0.4%

max_infected <- max(fit$I)
max_infected*0.004 
## [1] 58198.39

Peramalan Jumlah pada Kasus Positif Corona dengan Metode ARIMA

library(forecast)
## Warning: package 'forecast' was built under R version 3.5.3

Membuat data runtun waktu pada data.corona pada variabel Kasus.Kumulatif

data.corona <- ts(data.corona$Kasus..Kumulatif)

Visualisasi data.corona

autoplot(data.corona)

Membuat model ARIMA

model.arima <- auto.arima(data.corona, max.p = 5,max.q = 5,max.P = 2,max.Q = 2, max.d = 2, 
           max.D = 1,start.p = 0,start.q = 0,start.P = 0, start.Q = 0, 
           ic="aic", stepwise = F,approximation = F,test = "kpss",
           seasonal = F,lambda = "auto")

Melakukan peralamalan untuk 7 hari kedepan

forecast(model.arima,h=7)
##    Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 27       1183.314 1114.972 1253.871 1079.686 1292.122
## 28       1329.876 1209.963 1455.970 1148.970 1525.239
## 29       1485.781 1307.608 1676.385 1218.273 1782.361
## 30       1651.119 1406.163 1917.584 1285.097 2067.445
## 31       1825.977 1504.960 2180.968 1348.586 2382.824
## 32       2010.437 1603.618 2467.657 1408.332 2730.407
## 33       2204.581 1701.865 2778.663 1464.097 3111.952

Membuat visualisasi peramalan

autoplot(forecast(model.arima,h=7), ts.colour="orange")+ theme_gray() +xlab("Hari") +
  ylab("Jumlah Terinfeksi")