R/fitSIR.R

In cyhsuTN/COVID19: Prediction of Cases Volumne and Hospitalization by SIR Model

#' Estimation of cases volumne.
#' @description Cases volumne estimated from SIR model.
#' @import stats


#' @param susceptible Input. Number of susceptible people.
#' @param Infected Input. Number of infected people. If unknown,
#' it will be estimated by \emph{inihos}/((\emph{hosrate}/100) * (\emph{hms}/100)).
#' @param inihos Input. Initial number of hospitalized cases.
#' @param hosrate Input. Hospitalization rate of infected people (percentage between 0 to 100).
#' @param hms Input. Hospital market share (percentage between 0 to 100)
#' @param inidbt Input. Initial doubling time.
#' @param mrt Input. Mean recovery time.
#' @param sdis Input. Social distancing (percentage between 0 to 100).
#' @param dayFT Input. Days from today.

#' @return \item{result}{Cases volumne.}
#' @return \item{Mdbt}{Doubling time after mitigation.}
#' @return \item{EBRN}{Effective basic reproduction number.}


#' @examples ## To predicte 100 days from today (dayFT=100).
#' @examples casevolumne <- fitSIR(susceptible=4119405, Infected=3733, inihos=14,
#' @examples      hosrate=2.5, hms=15, inidbt=4, mrt=14, sdis=30, dayFT=100)
#' @examples head(casevolumne$result, 21) ## show the first 20 days
#'

#' @export


fitSIR <- function(susceptible, Infected, inihos=14, hosrate=2.5, hms=15, inidbt=4, mrt=14, sdis=30, dayFT=100) {

  S0 <- susceptible
  R0 <- 0
  if (missing(Infected)) I0 <- round(inihos/(hosrate/100)/(hms/100)) else I0 <- Infected

  douT <- inidbt
  sdis <- sdis

  gamma <- 1/mrt
  beta <- (2^(1/douT) - 1 + gamma)/S0
  betaS0 <- 2^(1/douT) - 1 + gamma

  # basic reproduction number, BRN = betaS0/gamma
  BRN <- betaS0/gamma
  EBRN <- BRN * (1-sdis/100) # effective basic reproduction number

  # EBRN*gamma = new betaS0;  new g = new betaS0 - gamma; to reestimate doubling time
  MdouT <- 1/log2(EBRN*gamma - gamma + 1)
  beta <- EBRN*gamma/S0

  SIRfun2 <- function(ts) {
    k <- 1
    reps <- matrix(NA, ts, 4)
    for (tk in 1:ts) {
      S1 <- (1 - beta * I0) * S0
      I1 <- (beta * S0 - gamma + 1) * I0
      R1 <- gamma * I0 + R0
      S0 <- S1
      I0 <- I1
      R0 <- R1
      reps[k,] <- round(c(k, S0, I0, R0))
      k <- k + 1
    }
    colnames(reps) <- c("days.from.today","susceptible","infected","recovered")
    reps
  }

  reports <- rbind(c(0, S0, I0, R0), SIRfun2(ts=dayFT))

  list(result=data.frame(reports), Mdbt=round(MdouT,2), EBRN=round(EBRN,2))

}