Pemodelan Kasus Persebaran COVID-19 Prov. Jawa Barat Dengan Metode SIR
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)
## Warning: package 'readxl' was built under R version 4.1.3
data.corona <- read_excel("covid19jabar.xlsx")
head(data.corona)
## # A tibble: 6 x 5
## tanggal kasus.kumulatif kasus.baru sembuh.baru meninggal.baru
## <dttm> <dbl> <dbl> <dbl> <dbl>
## 1 2020-12-01 00:00:00 1619 878 716 13
## 2 2020-12-02 00:00:00 2383 764 386 4
## 3 2020-12-03 00:00:00 4031 1648 253 21
## 4 2020-12-04 00:00:00 5023 992 379 4
## 5 2020-12-05 00:00:00 6109 1086 685 4
## 6 2020-12-06 00:00:00 7497 1388 1125 7
#sumber data: https://pikobar.jabarprov.go.id/data
#Membuat data harian dan menggabungkannya dalam data.corona
hari <- 1:31
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)
## Warning: package 'ggplot2' was built under R version 4.1.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 Jawa Barat")+theme(axis.text.x=element_text(angle = 45,hjust = 1))
plot.corona
#Membuat pemodelan penyebaran virus corona dengan model SIR berdasarkan jumlah penduduk Jawa Barat mengubah data ke tipe numeric.
Infected <- as.numeric(data.corona$kasus.kumulatif)
Day <- as.numeric(data.corona$hari)
N <- 49935858
R <- as.numeric(data.corona$sembuh.kumulatif)
#sumber data populasi: https://jabar.bps.go.id/
#Model SIR
library(deSolve)
## Warning: package 'deSolve' was built under R version 4.1.3
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.5612962 0.4387045
#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 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 Barat",outer=T,line=-2)
#Memanggil nilai R0
R0 <- setNames(opt_par["beta"]/opt_par["gamma"],"R0")
print(R0)
## R0
## 1.27944