R/hospfig.R

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

#' A figure for predicted hospitalization.
#' @description To show a figure for predicted hospitalization.
#' @import stats ggplot2
#' @importFrom reshape2 melt


#' @param showdays Input. Number of days to project.
#' @param hospitalization Input. Object from function \emph{Prehos}.
#' @param maxloadhos Input. Maximum load for hospitalized cases.
#' @param maxloadicu Input. Maximum load for ICU cases.
#' @param maxloadvent Input. Maximum load for ventilated cases.


#' @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 hospitalization <- Prehos(casevolumne, inihos=14, iniicu=0, iniven=0,
#' @examples      hosrate=2.5, icurate=0.75, venrate=0.5, outhosdays=7, outicudays=9,
#' @examples      outvendays=10, hms=15)
#' @examples hospfig(showdays=30, hospitalization,
#' @examples      maxloadhos=134, maxloadicu=72, maxloadvent=) ## show the first 30 days
#'

#' @export



hospfig <- function(showdays=60, hospitalization, maxloadhos=134, maxloadicu=72, maxloadvent) {
  if (showdays > nrow(hospitalization) - 1) {
    showdays <- nrow(hospitalization) - 1
  }

  dd <- reshape2::melt(hospitalization[1:(showdays+1),], id=c("days.from.today"))
  levels(dd$variable) <- c("Hospitalized", "ICU", "Ventilated")
  pp <- ggplot() +
    geom_line(data=dd, aes(x=days.from.today, y=value, colour=variable), size=1) +
    geom_point(data=dd, aes(x=days.from.today, y=value, colour=variable), size=2) +
    xlab("Date") +
    ylab("Patients in Hospita") +
    theme(axis.text=element_text(size=16), ## change size of text
          axis.title=element_text(size=18, face="bold")) + # change size of label name
    theme(legend.text=element_text(size=16)) + # change text size in legend
    theme(legend.title=element_text(size=16)) + # change title size of legend
    guides(col=guide_legend("Level")) +
    scale_colour_manual(values=c("blue","orange","red"))

  currentCDTDate <- function() as.Date(format(Sys.time(), tz="America/Chicago"))

  break1 <- round(seq(0,showdays,length.out=7))
  labels1 <- format(break1 + currentCDTDate(), "%b/%d")
  pp <- pp + scale_x_continuous(breaks=break1, labels=labels1) +
        theme(axis.text.x = element_text(size=14, angle=30, vjust = 0.8))


  if(!missing(maxloadhos)) {
    TF <- which(hospitalization[,"hosp"] > maxloadhos)
    if(sum(TF)==0) {
      daycrosshos <- Inf
    } else {
      daycrosshos <- min(TF) - 1
    }
    pp <- pp +
      geom_segment(aes(x = 0, y = maxloadhos, xend = min(daycrosshos, showdays), yend = maxloadhos), size=1, col="blue") +
      annotate(geom="text", x=0, y=maxloadhos+10, label=paste("Maxload =", maxloadhos),
               color="blue", size=5, hjust = 0, vjust = 0)
    if (daycrosshos <= showdays) {
      pp  <- pp + geom_segment(aes(x = daycrosshos, y = 0, xend = daycrosshos, yend = maxloadhos), size=1, col="blue") +
        annotate(geom="text",
                 x=daycrosshos, y=-1,
                 label=format(daycrosshos + currentCDTDate(), "%b/%d"),
                 color="blue", size=5)
    }
  }

  if(!missing(maxloadicu)) {
    TF <- which(hospitalization[,"icu"] > maxloadicu)
    if(sum(TF)==0) {
      daycrossicu <- Inf
    } else {
      daycrossicu <- min(TF) - 1
    }
    pp <- pp +
      geom_segment(aes(x = 0, y = maxloadicu, xend = min(daycrossicu, showdays), yend = maxloadicu), size=1, col="orange") +
      annotate(geom="text", x=0, y=maxloadicu+10, label=paste("Maxload =", maxloadicu),
               color="orange", size=5, hjust = 0, vjust = 0)
    if (daycrossicu <= showdays) {
      pp <- pp + geom_segment(aes(x = daycrossicu, y = 0, xend = daycrossicu, yend = maxloadicu), size=1, col="orange") +
        annotate(geom="text",
                 x=daycrossicu, y=-1,
                 label=format(daycrossicu + currentCDTDate(), "%b/%d"),
                 color="orange", size=5)
    }
  }

  if(!missing(maxloadvent)) {
    TF <- which(hospitalization[,"vent"] > maxloadvent)
    if(sum(TF)==0) {
      daycrossvent <- Inf
    } else {
      daycrossvent <- min(TF) - 1
    }
    pp <- pp +
      geom_segment(aes(x = 0, y = maxloadvent, xend = min(daycrossvent, showdays), yend = maxloadvent), size=1, col="red") +
      annotate(geom="text", x=0, y=maxloadvent+10, label=paste("Maxload =", maxloadvent),
               color="red", size=5, hjust = 0, vjust = 0)
    if (daycrossvent <= showdays) {
      pp <- pp + geom_segment(aes(x = daycrossvent, y = 0, xend = daycrossvent, yend = maxloadvent), size=1, col="red") +
        annotate(geom="text",
                 x=daycrossvent, y=maxloadvent-1,
                 label=format(daycrossvent + currentCDTDate(), "%b/%d"),
                 color="red", size=5)
    }
  }
  pp

}