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

UIN Maulana Malik Ibrahim Malang

#Model SIR adalah model matematika yang menggambarkan dinamika penyebaran penyakit menular (Ma dan Li, 2009:1) Pada model ini, populasi manusia dibagi menjadi tiga kelompok individu, yaitu kelompok individu rentan (Susceptible), kelompok individu terinfeksi (Infective), dan kelompok individu sembuh (Recovered).

#membaca data harian kasus covid-19 Jawa Tengah bulan Maret-Agustus 2020 dalam format excel

library(readxl)
data.corona <- read_excel("covidJateng.xlsx")

head(data.corona)
## # A tibble: 6 x 7
##   tanggal    hari positif meninggal sembuh positif_kumulatif sembuh_kumulatif
##   <chr>     <dbl>   <dbl>     <dbl>  <dbl>             <dbl>            <dbl>
## 1 44046         1       1         1      0                 1                0
## 2 44077         2       0         0      0                 1                0
## 3 44107         3       0         0      0                 1                0
## 4 44138         4       2         1      1                 3                1
## 5 44168         5       0         0      0                 3                1
## 6 3/13/2020     6       0         0      0                 3                1

#Visualisasi data kasus kumulatif covid-19 Jawa Tengah Maret-Agustus 2020

library(ggplot2)
hari <- data.corona$hari
y <- data.corona$positif_kumulatif
plot.corona <- ggplot(data = data.corona, aes(x=hari,y))+
  geom_point()+geom_smooth(method = "loess", formula = y~x,se=T)+ylab("Jumlah")+
  labs(title = "Total Infeksi COVID-19 Jawa Tengah Maret-Agustus 2020")+theme(axis.text.x=element_text(
    angle = 45,hjust = 1
  ))
plot.corona

#Membuat pemodelan penyebaran virus corona dengan model SIR berdasarkan jumlah penduduk Provinsi Jawa Tengah tahun 2020

Infected <- as.numeric(data.corona$positif_kumulatif)
Day <- as.numeric(data.corona$hari)
N <- 36516035
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*S*I/N
    dI <- beta*S*I/N-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.5310874 0.4689129

#Membuat data persebaran dengan waktu persebaran 360 hari

t <- 1:360
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 = "Hari Ke-", ylab="Jumlah Manusia",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 Tengah Maret-Agustus 2020",outer=T,line=-2)

#Memanggil nilai RO

fit$reff<- opt_par["beta"]/opt_par["gamma"]*fit$S/(fit$S+fit$I+fit$R)
fit$RO<- opt_par["beta"]/opt_par["gamma"]
RO <- setNames(opt_par["beta"]/opt_par["gamma"],"RO")
print(RO)
##       RO 
## 1.132593
ggplot() + geom_line (data=fit, aes(x=time, y= reff))+geom_point (data=fit, aes(x=time,y=reff), colour="red")+ geom_line (data=fit,aes(x=time, y=RO), colour="green")+xlab("Time (days)") + ylab ("Reff")+ labs (title = paste("Reproduction number levels: Beta =", opt_par["beta"],"and Gamma=", opt_par["gamma"]))