Dosen Pengampu : Prof. Dr. Suhartono, M.Kom.

UIN Maulana Malik Ibrahim Malang

Pemodelan Kasus Persebaran COVID-19 Dengan Model SIR di Jawa Timur

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 Kasus Harian

library(readxl)
data.corona <- read_excel("covid19jatim.xlsx")
## Warning in strptime(x, format, tz = tz): unable to identify current timezone 'C':
## please set environment variable 'TZ'
## New names:
## * `` -> ...6
## * `` -> ...7
## * `` -> ...8
head(data.corona)
## # A tibble: 6 x 8
##   tanggal             kasus.kumulatif kasus.baru sembuh.baru meninggal.baru
##   <dttm>                        <dbl>      <dbl>       <dbl>          <dbl>
## 1 2020-04-01 00:00:00             126          9          77              4
## 2 2020-04-02 00:00:00             126          0          91             12
## 3 2020-04-03 00:00:00             173         47           7              1
## 4 2020-04-04 00:00:00             178          5           3              2
## 5 2020-04-05 00:00:00             213         35           8              2
## 6 2020-04-06 00:00:00             222          9          11              5
## # ... with 3 more variables: ...6 <lgl>, ...7 <dttm>, ...8 <dbl>

sumber data: https://www.kaggle.com/datasets/hendratno/covid19-indonesia

Membuat data harian dan menggabungkannya dalam data.corona

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

Membuat data kumulatif pada kolom sembuh baru

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

Visualisasi kasus.kumulatif

library(ggplot2)
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 Jawa Timur")+theme(axis.text.x=element_text(angle = 45,hjust = 1))
plot.corona

Membuat pemodelan penyebaran virus corona dengan model SIR berdasarkan jumlah penduduk Jawa Timur mengubah data ke tipe numeric.

Infected <- as.numeric(data.corona$kasus.kumulatif)
Day <- as.numeric(data.corona$hari)
N <- 10609681
R <- as.numeric(data.corona$sembuh.kumulatif)

sumber data populasi: https://jawatimur.bps.go.id/

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] "ERROR: ABNORMAL_TERMINATION_IN_LNSRCH"

Memanggil nilai beta dan gamma

opt_par <- setNames(opt$par,c("beta","gamma"))
opt_par
##      beta     gamma 
## 0.5307859 0.4744822

Membuat data persebaran dengan waktu persebaran 150 hari

t <- 1:240
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 Orang",lwd = 2,lty=1,col=col,log="y")
## Warning in xy.coords(x, y, xlabel, ylabel, log = log, recycle = TRUE): 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 Jawa Timur",outer=T,line=-2)

Memanggil nilai R0

R0 <- setNames(opt_par["beta"]/opt_par["gamma"],"R0")
print(R0)
##       R0 
## 1.118663

Peramalan Jumlah pada Kasus Positif Corona dengan Metode ARIMA

library(forecast)
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo

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 30 hari kedepan

forecast(model.arima,h=30)
##     Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
##  92       12100.48 11873.33 12330.22 11754.13 12452.88
##  93       12445.92 12126.84 12769.98 11959.93 12943.55
##  94       12758.99 12314.14 13213.32 12082.47 13457.68
##  95       13076.72 12521.93 13645.97 12234.05 13953.20
##  96       13399.14 12741.24 14076.96 12400.97 14443.91
##  97       13726.28 12968.41 14510.04 12577.61 14935.52
##  98       14058.14 13201.55 14947.16 12761.08 15431.03
##  99       14394.76 13439.51 15389.52 12949.65 15932.29
## 100       14736.16 13681.59 15837.93 13142.20 16440.51
## 101       15082.36 13927.26 16292.95 13337.96 16956.56
## 102       15433.38 14176.18 16755.03 13536.39 17481.11
## 103       15789.24 14428.07 17224.48 13737.06 18014.68
## 104       16149.97 14682.73 17701.61 13939.65 18557.69
## 105       16515.59 14939.99 18186.64 14143.91 19110.50
## 106       16886.12 15199.71 18679.77 14349.65 19673.43
## 107       17261.57 15461.79 19181.20 14556.68 20246.76
## 108       17641.98 15726.13 19691.07 14764.87 20830.73
## 109       18027.36 15992.66 20209.54 14974.10 21425.57
## 110       18417.74 16261.32 20736.75 15184.26 22031.50
## 111       18813.13 16532.03 21272.83 15395.27 22648.72
## 112       19213.55 16804.75 21817.90 15607.04 23277.41
## 113       19619.04 17079.44 22372.08 15819.51 23917.76
## 114       20029.60 17356.05 22935.47 16032.61 24569.95
## 115       20445.26 17634.55 23508.20 16246.29 25234.13
## 116       20866.04 17914.90 24090.36 16460.50 25910.48
## 117       21291.96 18197.07 24682.06 16675.18 26599.16
## 118       21723.04 18481.03 25283.40 16890.29 27300.32
## 119       22159.30 18766.76 25894.48 17105.80 28014.11
## 120       22600.77 19054.22 26515.39 17321.66 28740.70
## 121       23047.45 19343.40 27146.24 17537.84 29480.23

Membuat visualisasi peramalan

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