Author : Lia Wahyuliningtyas

library(readxl)
data.corona <- read_excel(path = "C:/Users/Windows 10/Desktop/kotabanjarmasin2.xlsx")
data.corona
## # A tibble: 55 x 4
##    tanggal             Positif Sembuh Meninggal
##    <dttm>                <dbl>  <dbl>     <dbl>
##  1 2020-08-13 00:00:00    2542   1311       146
##  2 2020-08-14 00:00:00    2555   1311       147
##  3 2020-08-15 00:00:00    2573   1330       147
##  4 2020-08-16 00:00:00    2581   1330       147
##  5 2020-08-17 00:00:00    2609   1411       147
##  6 2020-08-18 00:00:00    2625   1411       147
##  7 2020-08-19 00:00:00    2657   1482       147
##  8 2020-08-20 00:00:00    2687   1535       149
##  9 2020-08-21 00:00:00    2690   1804       149
## 10 2020-08-22 00:00:00    2690   1804       149
## # ... with 45 more rows
hari <- 1:55
data.corona <- cbind(data.corona,hari)

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

library(ggplot2)

plot.corona <- ggplot(data = data.corona, aes(x=hari,y=Positif))+
  geom_point()+geom_smooth(method = "loess", formula = y~x,se=T)+ylab("Jumlah")+
  labs(title = "Total Infeksi COVID-19 Kota Banjarmasin")+theme(axis.text.x=element_text(
    angle = 45,hjust = 1
  ))
plot.corona

Infected <- as.numeric(data.corona$Positif)
Day <- as.numeric(data.corona$hari)
N <- 667000
R <-as.numeric(data.corona$Sembuh)

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"
opt_par <- setNames(opt$par,c("beta","gamma"))
opt_par
##      beta     gamma 
## 0.5145690 0.4854998
t <- 1:90
fit <- data.frame(ode(y=init, times=t,func=SIR, parms = opt_par))

col <- 1:4
matplot(fit$time, fit[,2:4],type="l",xlab = "Hari ke-", ylab="Jumlah Penduduk",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 Kota Banjarmasin Kalsel",outer=T,line=-2)

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