R/covid_simr-gamma.R

In nhs-bnssg-analytics/covid-simr: easy-to-use function for modelling the effect of different mitigating scenarios on projected capacity-dependent covid19 deaths

#' Covid-19 ICU bed-occupancy simulation
#'
#' @param cases \code{data.frame}, hosptialisation data, one column named \code{date} (string), and one column named \code{hospitalisations} which is an integer recording new cases for that date.
#' @param los_median \code{positive numeric}, median length of stay (in days) for admitted covid19 pateints, to be estimated by user from data or other source
#' @param los_95 \code{positive numeric}, 0.95 quantile length of stay (in days) for admitted covid19 patients, to be estimated by user from data or other source
#' @param cap \code{positive integer}, maximum number of patients who can be concurrently admitted to the hospital/unit
#' @param pfat \code{numeric} in range 0 (no chance of death) to 1 (certain death), probability that a patient who arrives but is rejected/cannot be admitted (because there are no beds free) dies
#' @param tol \code{numeric} between 0 and 100, with 0 corresponding to full confidence inputs are correct, with 100 least level of confidence, subjective assessment of reliability of input hospitalisation estiamtes
#' @param nreps \code{positive integer}, number of simulation replications to perform - larger number means better results, but longer time to compute
#'
#' @return An object of type \code{Covidsimr}, which is a list including the slots:
#'   \itemize{
#'     \item{\code{data},}{ Output results from simulation}
#'     \item{\code{data_cum},}{ Output results from simulation, in cumulative form}
#'     \item{\code{tol},}{ \code{TRUE/FALSE} based on whether a tolerance has been set in the function inputs}
#'     \item{\code{cap},}{ Maximum ICU capactiy}
#'   }
#'   Of most interest to users are \code{data} and \code{data_cum}, the other elements are included for use in the plotting functions.
#' @export
#'
#' @importFrom magrittr %>%
#'
#' @examples
covid_simr <- function(cases,
                       los_median,
                       los_95,
                       cap,
                       pfat,
                       tol = 25,
                       nreps = 100) {

  # pre-alloacte symbols to NULL to avoid NOTES from CMD check
  dates <- metric <- value <- fill <- time <- outcome <- type <- value_cum <-
    q01 <- q025 <- q05 <- q15 <- q30 <- q70 <- q85 <- q95 <- q975 <- q99 <-
    hospitalisations <- occ_q01 <- occ_q025 <- occ_q05 <- occ_q15 <- occ_q30 <-
    occ_q70 <- occ_q85 <- occ_q95 <- occ_q975 <- occ_q99 <- occ_median <-
    rejected_died_q01 <- rejected_died_q025 <- rejected_died_q05 <-
    rejected_died_q15 <- rejected_died_q30 <- rejected_died_q70 <-
    rejected_died_q85 <- rejected_died_q95 <- rejected_died_q975 <-
    rejected_died_q99 <- rejected_died_median <-
    admitted_q01 <- admitted_q025 <- admitted_q05 <- admitted_q15 <-
    admitted_q30 <- admitted_q70 <- admitted_q85 <- admitted_q95 <-
    admitted_q975 <- admitted_q99 <- admitted_median <-
    rejected_survived_q01 <- rejected_survived_q025 <- rejected_survived_q05 <-
    rejected_survived_q15 <- rejected_survived_q30 <- rejected_survived_q70 <-
    rejected_survived_q85 <- rejected_survived_q95 <- rejected_survived_q975 <-
    rejected_survived_q99 <- rejected_survived_median <- NULL


  # calculate lognormal distribution parameters from given information (los_mean and los_95)
  meanlog<-log(los_median)
  sdlog<-((sqrt(2)*stats::qnorm((1+0.9)/2)/sqrt(2))^-1)*log(los_95/los_median)

  # determine whether epidemic curve uncertainty is to be considered
  if (tol > 0 & tol <= 100) {
    uncert <- TRUE
  } else {
    uncert <- FALSE
  }


  #This section Reads in given daily hospitalisations esimtate from a CSV (based on epidemic curve data, not caculated by this package).
  #    Generates confidence intervals by sampling noise from normal distrubtion, scaling that by a tolerance parameter (given as
  #    a percentage uncertainly as an integer in the range 0-100)  and scaling the orignial hospitalisations estimate using
  #    that (i.e. orignial value * 1+scaling, where the scaling can be +-, floored at -1). The output of this function is
  #    then used internally as the input dataframe by the run_sim_fun function in this package.


  #create three column dataframe:
  #date (date format)
  #name of metric type (e.g. hospitalisations) (character format)
  #value for date/type pair (numeric format)
  inputs <- cases
  inputs <- inputs %>%
    dplyr::mutate(dates = as.Date(dates, "%d/%m/%Y")) %>%
    tidyr::pivot_longer(cols = -dates,
                        names_to = "metric",
                        values_to = "value")

  #If a nonzero tolerance has been set (strictly between 0 and 100 - for percent)
  #then uncertaintly about the epidemic curve is assumed
  #and daily arrivals are modified by adding some normally distrubuted noise
  #(floored at 0, and increased by a maximum of the upper positive limit of a sample from
  #the normal distrubtion about zero with sd 1/3), modified by choice of tolerance parameter (to sim, not normal)
  #For each date, is is repeated reps number of times and the results use to derive quantile estimates for possible arrivals
  #to be used as confidence intervals on the original central estimate
  if (uncert == TRUE) {
    inputs_sim <- lapply(1:nreps, function(i) {
      set.seed(i)
      scaler <-
        tol / 100 * max(stats::rnorm(1, mean = 0, sd = 1 / 3),-1)
      data.frame(
        dates = inputs$dates,
        metric = paste0("rep", i),
        value = inputs$value * (1 + scaler)
      )
    })
    inputs_sim <- do.call("rbind", inputs_sim)
    inputs <- rbind(inputs, inputs_sim)
    #get CIs
    inputs_sim_ci <- inputs %>%
      dplyr::group_by(dates) %>%
      dplyr::filter(substr(metric, 1, 3) == "rep") %>%
      dplyr::summarise(
        q30 = stats::quantile(value, 0.3),
        q15 = stats::quantile(value, 0.15),
        q05 = stats::quantile(value, 0.05),
        q025 = stats::quantile(value, 0.025),
        q01 = stats::quantile(value, 0.01),
        q70 = stats::quantile(value, 0.7),
        q85 = stats::quantile(value, 0.85),
        q95 = stats::quantile(value, 0.95),
        q975 = stats::quantile(value, 0.975),
        q99 = stats::quantile(value, 0.99)
      ) %>%
      tidyr::pivot_longer(cols = -dates,
                          names_to = "metric",
                          values_to = "value")
    #adds original inputs datframe now has multiple rows for each date
    #one of orignial hospitalisations estimate,
    #for each rep, one with a modified hospitalisations estimate (scaled by the normal sample - itself scaled by the tolerance value)
    #and one for each of the confidence quantiles derived rmo the modified hospitalisations
    inputs <- rbind(inputs, inputs_sim_ci)
  }

  #Now run nreps replications of the simulation function to get results for nreps possible
  #actualisations of the parameters

  #Do the runs in parallel to speed up processing time
  #Set up the clusters (set to one fewer than the total number of logical cores in the user's comput)
  cl <- parallel::makeCluster(parallel::detectCores() - 1)
  #send necessary objects created in the global environment to the clusters
  parallel::clusterExport(
    cl = cl,
    varlist = c("uncert", "inputs", "los_median", "los_95", "cap", "pfat"),
    envir = environment()
  )
  #make the libraries used within the simfn function available to the clusters
  parallel::clusterEvalQ(cl = cl,
                         c(library(tidyr),
                           library(dplyr)))
  #return results of the replications
  RES <- parallel::parLapply(cl, 1:nreps,
                             simfn,
                             dates = inputs$dates,
                             metric = inputs$metric,
                             value = inputs$value,
                             los_median = los_median,
                             los_95 = los_95,
                             cap = cap,
                             pfat = pfat,
                             uncert = uncert)
  parallel::stopCluster(cl)

  #combine the simulation outputs into a single dataframe
  outputs_sim <- do.call("rbind", RES)
  outputs <- outputs_sim %>%
    tidyr::pivot_longer(
      cols = -c(rep, time),
      names_to = "outcome",
      values_to = "value"
    ) %>%
    dplyr::group_by(time, outcome) %>%
    dplyr::summarise(
      mean = mean(value),
      median = stats::median(value),
      q30 = stats::quantile(value, 0.3),
      q25 = stats::quantile(value, 0.25),
      q15 = stats::quantile(value, 0.15),
      q05 = stats::quantile(value, 0.05),
      q025 = stats::quantile(value, 0.025),
      q01 = stats::quantile(value, 0.01),
      q70 = stats::quantile(value, 0.7),
      q75 = stats::quantile(value, 0.75),
      q85 = stats::quantile(value, 0.85),
      q95 = stats::quantile(value, 0.95),
      q975 = stats::quantile(value, 0.975),
      q99 = stats::quantile(value, 0.99)
    ) %>%
    tidyr::pivot_longer(
      cols = -c(time, outcome),
      names_to = "type",
      values_to = "value"
    ) %>%
    tidyr::pivot_wider(names_from = c(outcome, type),
                       values_from = value) %>%
    dplyr::mutate(dates = seq(
      from = min(inputs$dates) + time - 1,
      by = 1,
      length.out = 1
    ))

  data_cum <- outputs_sim %>%
    tidyr::pivot_longer(
      cols = -c(rep, time),
      names_to = "outcome",
      values_to = "value"
    ) %>%
    dplyr::group_by(outcome, rep) %>%
    dplyr::mutate(value_cum = cumsum(value)) %>%
    dplyr::select(-value) %>%
    dplyr::group_by(time, outcome) %>%
    dplyr::summarise(
      mean = mean(value_cum),
      median = stats::median(value_cum),
      q30 = stats::quantile(value_cum, 0.3),
      q25 = stats::quantile(value_cum, 0.25),
      q15 = stats::quantile(value_cum, 0.15),
      q05 = stats::quantile(value_cum, 0.05),
      q025 = stats::quantile(value_cum, 0.025),
      q01 = stats::quantile(value_cum, 0.01),
      q70 = stats::quantile(value_cum, 0.7),
      q75 = stats::quantile(value_cum, 0.75),
      q85 = stats::quantile(value_cum, 0.85),
      q95 = stats::quantile(value_cum, 0.95),
      q975 = stats::quantile(value_cum, 0.975),
      q99 = stats::quantile(value_cum, 0.99)
    ) %>%
    tidyr::pivot_longer(
      cols = -c(time, outcome),
      names_to = "type",
      values_to = "value_cum"
    ) %>%
    tidyr::pivot_wider(names_from = c(outcome, type),
                       values_from = value_cum) %>%
    dplyr::mutate(dates = seq(
      from = min(inputs$dates) + time - 1,
      by = 1,
      length.out = 1
    ))


  ####################################################################################################################################################################

  #save inputs (with tolerances/confinece range) and outputs as csvs,
  #in same location script was run from
  #write.csv(inputs,paste0("inputs",filename_ext,".csv"),row.names=FALSE)
  # utils::write.csv(
  #   cbind(data.frame(cap = cap), as.data.frame(outputs)),
  #   paste0("outputs", filename_ext, ".csv"),
  #   row.names = FALSE
  # )
  # utils::write.csv(
  #   cbind(data.frame(cap = cap), as.data.frame(data_cum)),
  #   paste0("data_cum", filename_ext, ".csv"),
  #   row.names = FALSE
  # )

  #return outputs
  out <- structure(list(data = outputs,
                        data_cum = data_cum,
                        inputs = inputs,
                        uncert = uncert,
                        cap = cap),
                   class = "Covidsimr")
  return(out)
}

#' COVID-19 ICU bed occupancy simulation
#'
#' @param rep (\code{postive integer}), number of repetitions in simulation.
#' @param dates (\code{vec, Date}), dates in siulation.
#' @param metric PLACEHOLDER
#' @param value PLACEHOLDER
#' @param uncert \code{log}, does the simulation include uncertainty.
#' @inheritParams covid_simr
#'
#' @return Results of simulation.
#'
#' @examples
simfn <- function(rep, dates, metric, value, los_median, los_95, cap, pfat, uncert) {
  # calculate lognormal distribution parameters from given information (los_mean and los_95)
  meanlog<-log(los_median)
  sdlog<-((sqrt(2)*qnorm((1+0.9)/2)/sqrt(2))^-1)*log(los_95/los_median)
  #fix the random number stream for this replication
  #(so that if run again without changing any substantive parameters, the results will be the same)
  #uses rep number so that a different random number stream is initialised for each replication
  set.seed(rep)
  dmax <- length(unique(dates))
  #get num arrivals by day
  #then spread them out over time within days
  arr_times <-
    unlist(sapply(1:dmax, function(x) {
      #hospitalisation estimate for each day
      #either the original value from the epidemic curve
      #or, if a tolrenace value has been set, the modified value for the given replication
      narr <-
        value[which(metric == ifelse(uncert == TRUE, paste0("rep", rep), "hospitalisations"))][x]
      #For each date, a (re)sampled total hospitalisation number, using the original value as the rate parameter
      #to the Poisson distribution
      narr.pois <- stats::rpois(1, lambda = narr)
      #spread these arrivals out randomly throughout the given day, ordering first to last
      sort(stats::runif(narr.pois, 0, 1) + x - 1)
    }))
  #"cal" - the simulation schedule
  #used internally within the simulation loop to track events and determine admissions, rejections, service ends
  #dataframe listing the arrivals for the given day in order:
  #"id" - positive integer, index column
  #"time" - numeric, decimal time between first and last arrivals
  #"event" - character, describes type of event - intinally, here, that is just a patient arrival
  cal <-
    data.frame(id = 1:length(arr_times),
               time = arr_times,
               event = "arrival")
  #The simulation clock - used to update both the simulation schedule "cal" and (by rounding up) the outputs dataframe "res"
  tx <- 0
  #initialise outcome dataframe "res" for each day
  #number of rows is approx number of days plus 1/3
  #each row tracks the state of the hospital in terms of total occupancy at, and total admissions and rejections (separated into
  #survived and died) by that time index
  #dataframe cols:
  #"time" integer, tracks days from start to end + ~1/3
  #"occ", integer, current total numeber of beds occupied
  #"admitted", integer, cumulative number patients admitted since start
  #"rejected_died", integer, cumulative number of patients who were refused admission and as a result died, since start
  #"rejected_survived", integer, cumulative number of patients who were refused admission but survived, since start
  #CHANGE - initialise results with occupancy NA instead of occupancy zero - then handle occupancies at end ####
  res <-
    data.frame(
      time = 1:round(dmax * 1.33),
      occ = NA,
      admitted = 0,
      rejected_died = 0,
      rejected_survived = 0
    )
  #initialise hospital occupancy before starting simulation loop ####
  occ <- 0  #occupancy (number in unit)
  while (nrow(cal) > 0 &
         tx < max(res$time)) {
    #find all rows in the schedule which have not happened yet (time greater than simulation clock time tx)
    #and which are either arrivals or service completions
    #return as a vector of row indices to the schedule dataframe "cal"
    ind1 <- which(cal$time > tx &
                    cal$event %in% c("arrival", "endsrv"))
    #find the row of the schedule which corresponds to the earliest of those events and return its index
    ind <- ind1[which.min(cal$time[ind1])]
    tx_old<-tx
    occ_old<-occ
    #advance the simulation clock to that time
    tx <- cal$time[ind]
    #return the start time of the next full schedule day after the current simulation clock time
    tx_day <- ceiling(tx)
    #if the next day is beyond the time limit (the maximum arrivals schedule date plus approx 1/3), then stop the simulation loop
    #,and proceed to return the results dataframe as the function output
    if (tx_day > max(res$time))
      break
    #otherwise, process events 3 types (1.arrival-admitted, 2.arrival-rejected[died/survived], and 3.end of service)
    #Arrival
    if (cal$event[ind] == "arrival") {
      #1. Arrival-admission: if there is a free bed, they are admitted to the hospital
      if (occ < cap) {
        #increment the cumulative tally of admitted patients at the start of the next full calendar date by one in the results
        #dataframe. This is an output measure and will not be used to decide whether or not subsequent arrivals can be admitted
        res$admitted[tx_day] <-
          res$admitted[tx_day] + 1
        #add the start of this patient's length of stay as an event, to the end of the simulation schedule, at the decimal start time
        #(n.b. not the rounded up to start of next day time used to track occupancy)
        #POTENTIALLY SLOW, AND DONE LOTS OF TIMES ####
        cal <-
          rbind(cal,
                data.frame(
                  id = cal$id[ind],
                  time = tx,
                  event = "startsrv"
                ))
        #for the patient just admitted, sample a length of stay from the lognormal distribution using the function input parameters,
        #which scale it to expected median length of stay for all patients (dervived externally by user
        #from case data or other source)
        los <-
          stats::rlnorm(1, meanlog = meanlog, sdlog = sdlog)
        #add the end of this patient's length of stay as an event to the end of the simulation schedule, specifying the event type
        #and the decimal time (given by the patient's service start time plus the LOS just sampled)
        #POTENTIALLY SLOW, AND DONE LOTS OF TIMES ####
        cal <-
          rbind(cal,
                data.frame(
                  id = cal$id[ind],
                  time = tx + los,
                  event = "endsrv"
                ))
        #increment the current actual occupancy (at this point in time by the simulation clock tx) by one
        #this will be used to check whether the next arrival can be admitted
        occ <- occ + 1
      } else {
        ##2. Arrival-rejection: There is not a free bed, so the patient is rejected
        #remove the event row from the simulation schedule
        #POTENTIALLY SLOW, AND DONE LOTS OF TIMES ####
        cal <-
          cal[-which(cal$id == cal$id[ind]),]
        #Does the rejected patient die? Sample from the range [0,1] with each value being equally probable. If the sample is less
        #than the user-inputed estiamte of a rejected patient dying, they die
        if (stats::runif(1, 0, 1) < pfat) {
          #2.1 Arrival-rejection-died
          #increment the running tally of patients who have been rejected and died (at the start of the next full day, in
          #the results dataframe). This is an output, it is not used to decide future events within the simulation loop.
          res$rejected_died[tx_day] <-
            res$rejected_died[tx_day] + 1
        } else {
          #2.2 Arrival-rejection-survived
          #sample from unif was greater than the given chance of dying, so patient survives
          #increment the running tally fo rejected-survived in the outputs dataframe (res) a the row corresponding to the start of
          #the next day. Thisis an output, it is not used to decide future events within the simulation loop.
          res$rejected_survived[tx_day] <-
            res$rejected_survived[tx_day] + 1
        }
      }
      #3. End of service
      #Note that this may be survival or death - but not death due to capacity constraint
      #i.e. it is normal outcome in event of sufficient capacity for this patient
    } else if (cal$event[ind] == "endsrv") {
      #Remove row for this patient from the simulation schedule
      #POTENTIALLY SLOW, AND DONE LOTS OF TIMES ####
      cal <-
        cal[-which(cal$id == cal$id[ind]),]    # this is no longer needed since p measures recorded separately at each iteration
      #de-increment occuapncy at current simulation time
      occ <- occ - 1
    }
    #sort the simulation schedule into ascending time order, by the time column
    #THIS WILL BE SLOW - AND DONE LOTS OF TIMES ####
    cal <- cal[order(cal$time),]
    #save results, extract performance measures
    res_ind<-which(res$time==tx_day)
    wt_new<-(tx-tx_old)/tx
    res$occ[res_ind]<-ifelse(is.na(res$occ[res_ind]),(tx-floor(tx))*occ_old+(ceiling(tx)-tx)*occ,wt_new*occ+(1-wt_new)*res$occ[res_ind])
  }
  #for time periods where there was no event, the occupancy field in the results will be NA
  #replace this to be (i) zero if it is the first day
  #                   (i) the same occupancy as the last non-NA day otherwise
  if (is.na(res$occ[1])) {
    res$occ[1] <- 0
  }
  res <- res %>% tidyr::fill(occ)
  return(cbind(data.frame(rep = rep), res))
}

#' Tidy print \code{Covidsimr} result.
#'
#' @param x a \code{covidsimR} object.
#' @param ... any other parameters
#'
#' @return a tidy print of the return data.
#' @export
#'
#' @examples
#' @rdname print
#' @method print Covidsimr
print.Covidsimr <- function(x, ...){
  print(tibble::tibble(x$data))
}

#' Plot output of COVID-19 simulation
#'
#' @param x a \code{CovidsimR} object.
#' @param ... any other parameters.
#' @export
#'
#' @return Mosaic plot of diagnostic charts from simulation.
#'
#'
#' @examples
plot.Covidsimr <- function(x, ...) {
  inputs  <- x$inputs
  outputs <- x$data
  data_cum <- x$data_cum
  uncert <- x$uncert
  cap <-x$cap

  #plot the results

  #daily hospitalisations (including resampling and tolerance to give confidence ranges)
  plot1 <- plot_daily_hospitalisations(inputs, outputs, uncert, display = FALSE)
  #bed occupancy over simulation period
  plot2 <- plot_daily_bed_occupancy(outputs, cap, display = FALSE)
  #deaths resulting from insufficient capacity over simulation period
  plot3 <- plot_daily_deaths(outputs, display = FALSE)
  #cumulative total of admitted patients over simulation period
  plot4 <- plot_cum_admissions(data_cum, display = FALSE)
  #cumulative total of patients who could not be admitted because of capacity constraints but did
  #surivive, over simulation period
  plot5 <- plot_cum_total_nonadmitted_surv(data_cum, display = FALSE)
  #cumulative total of patients who died as a result of insufficent capacity over simulation period
  plot6 <- plot_cum_died_inscap(data_cum, display = FALSE)
  gridExtra::grid.arrange(plot1, plot2, plot3, plot4, plot5, plot6, nrow = 2)

}




#' Plot daily hospitalisations (including resampling and tolerance to give confidence ranges)
#'
#' @inheritParams simfn
#' @param inputs  PLACEHOLDER
#' @param outputs \code{data.frame}, output data from simulation.
#' @param display \code{log}, defaults to \code{TRUE}. Should the plot display or not.
#'
#' @return Time series plot of daily hospitalistions. Invisibly returns the plot object.
#' @export
#'
#' @examples
plot_daily_hospitalisations <- function(inputs, outputs, uncert, display = TRUE) {
  #daily hospitalisations (including resampling and tolerance to give confidence ranges)
  plot <- inputs %>%
    tidyr::pivot_wider(names_from = "metric", values_from = "value") %>%
    ggplot2::ggplot(ggplot2::aes_string(x = "dates")) +
    ggplot2::labs(title = "Inputs: Cases requiring hospitalisation (per day)") +
    ggplot2::theme(
      axis.title.x = ggplot2::element_blank(),
      axis.title.y = ggplot2::element_blank(),
      plot.title = ggplot2::element_text(size = 11),
      axis.text.x = ggplot2::element_text(angle = 45, hjust = 1)
    ) +
    ggplot2::scale_x_date(
      date_breaks = "months",
      date_labels = "%b-%y",
      limits = c(min(outputs$dates), max(outputs$dates))
    )
  if (uncert == TRUE) {
    plot <- plot +
      ggplot2::geom_ribbon(ggplot2::aes_string(ymin = "q01", ymax = "q025"),
                           fill = "grey",
                           alpha = 0.6) +
      ggplot2::geom_ribbon(ggplot2::aes_string(ymin = "q025", ymax = "q05"),
                           fill = "dodgerblue1",
                           alpha = 0.6) +
      ggplot2::geom_ribbon(ggplot2::aes_string(ymin = "q05", ymax = "q15"),
                           fill = "dodgerblue2",
                           alpha = 0.6) +
      ggplot2::geom_ribbon(ggplot2::aes_string(ymin = "q15", ymax = "q30"),
                           fill = "dodgerblue3",
                           alpha = 0.6) +
      ggplot2::geom_ribbon(ggplot2::aes_string(ymin = "q30", ymax = "q70"),
                           fill = "dodgerblue4",
                           alpha = 0.6) +
      ggplot2::geom_ribbon(ggplot2::aes_string(ymin = "q70", ymax = "q85"),
                           fill = "dodgerblue3",
                           alpha = 0.6) +
      ggplot2::geom_ribbon(ggplot2::aes_string(ymin = "q85", ymax = "q95"),
                           fill = "dodgerblue2",
                           alpha = 0.6) +
      ggplot2::geom_ribbon(ggplot2::aes_string(ymin = "q95", ymax = "q975"),
                           fill = "dodgerblue1",
                           alpha = 0.6) +
      ggplot2::geom_ribbon(ggplot2::aes_string(ymin = "q975", ymax = "q99"),
                           fill = "grey",
                           alpha = 0.6)
  }
  plot <- plot +
    ggplot2::geom_line(ggplot2::aes_string(y = "hospitalisations"))

  if(display){print(plot)}
  invisible(plot)
}


#' Plot daily bed occupancy over simulation period
#'
#' @inheritParams plot_daily_hospitalisations
#' @param outputs
#' @param cap \code{int}, the bed capacity.
#'
#' @return Time-siries plot of daily bed occupancy. Invisibly returns the plot object.
#' @export
#'
#' @examples
plot_daily_bed_occupancy <- function(outputs, cap, display=TRUE) {
  #bed occupancy over simulation period
  plot <- outputs %>%
    ggplot2::ggplot(ggplot2::aes_string(x = "dates")) +
    ggplot2::geom_ribbon(ggplot2::aes_string(ymin = "occ_q01", ymax = "occ_q025"),
                         fill = "darkgray",
                         alpha = 0.6) +
    ggplot2::geom_ribbon(ggplot2::aes_string(ymin = "occ_q025", ymax = "occ_q05"),
                         fill = "darkorchid1",
                         alpha = 0.6) +
    ggplot2::geom_ribbon(ggplot2::aes_string(ymin = "occ_q05", ymax = "occ_q15"),
                         fill = "darkorchid2",
                         alpha = 0.6) +
    ggplot2::geom_ribbon(ggplot2::aes_string(ymin = "occ_q15", ymax = "occ_q30"),
                         fill = "darkorchid3",
                         alpha = 0.6) +
    ggplot2::geom_ribbon(ggplot2::aes_string(ymin = "occ_q30", ymax = "occ_q70"),
                         fill = "darkorchid4",
                         alpha = 0.6) +
    ggplot2::geom_ribbon(ggplot2::aes_string(ymin = "occ_q70", ymax = "occ_q85"),
                         fill = "darkorchid3",
                         alpha = 0.6) +
    ggplot2::geom_ribbon(ggplot2::aes_string(ymin = "occ_q85", ymax = "occ_q95"),
                         fill = "darkorchid2",
                         alpha = 0.6) +
    ggplot2::geom_ribbon(ggplot2::aes_string(ymin = "occ_q95", ymax = "occ_q975"),
                         fill = "darkorchid1",
                         alpha = 0.6) +
    ggplot2::geom_ribbon(ggplot2::aes_string(ymin = "occ_q975", ymax = "occ_q99"),
                         fill = "darkgray",
                         alpha = 0.6) +
    ggplot2::geom_hline(yintercept = cap,
                        linetype = "dashed",
                        colour = "darkgray") +
    ggplot2::geom_line(ggplot2::aes_string(y = "occ_median")) +
    ggplot2::labs(title = paste0("Outputs: Occupied beds (capacity inputted=", cap, ")")) +
    ggplot2::theme(
      axis.title.x = ggplot2::element_blank(),
      axis.title.y = ggplot2::element_blank(),
      plot.title = ggplot2::element_text(size = 11),
      axis.text.x = ggplot2::element_text(angle = 45, hjust = 1)
    ) +
    ggplot2::scale_x_date(date_breaks = "months", date_labels = "%b-%y")

  if(display){print(plot)}
  invisible(plot)
}

#' Plot daily deaths result from insufficient capacity over simulation period
#'
#' @inheritParams plot_daily_hospitalisations
#'
#' @return Time-series plot of daily deaths from insufficient capacity. Invisibly returns the plot object.
#' @export
#'
#' @examples
plot_daily_deaths <- function(outputs, display = TRUE) {
  #deaths resulting from insufficient capacity over simulation period
  plot <- outputs %>%
    ggplot2::ggplot(ggplot2::aes_string(x = "dates")) +
    ggplot2::geom_ribbon(
      ggplot2::aes_string(ymin = "rejected_died_q01", ymax = "rejected_died_q025"),
      fill = "yellow1",
      alpha = 0.6
    ) +
    ggplot2::geom_ribbon(
      ggplot2::aes_string(ymin = "rejected_died_q025", ymax = "rejected_died_q05"),
      fill = "orange2",
      alpha = 0.6
    ) +
    ggplot2::geom_ribbon(
      ggplot2::aes_string(ymin = "rejected_died_q05", ymax = "rejected_died_q15"),
      fill = "orangered3",
      alpha = 0.6
    ) +
    ggplot2::geom_ribbon(
      ggplot2::aes_string(ymin = "rejected_died_q15", ymax = "rejected_died_q30"),
      fill = "red3",
      alpha = 0.6
    ) +
    ggplot2::geom_ribbon(
      ggplot2::aes_string(ymin = "rejected_died_q30", ymax = "rejected_died_q70"),
      fill = "red2",
      alpha = 0.6
    ) +
    ggplot2::geom_ribbon(
      ggplot2::aes_string(ymin = "rejected_died_q70", ymax = "rejected_died_q85"),
      fill = "red3",
      alpha = 0.6
    ) +
    ggplot2::geom_ribbon(
      ggplot2::aes_string(ymin = "rejected_died_q85", ymax = "rejected_died_q95"),
      fill = "orangered3",
      alpha = 0.6
    ) +
    ggplot2::geom_ribbon(
      ggplot2::aes_string(ymin = "rejected_died_q95", ymax = "rejected_died_q975"),
      fill = "orange2",
      alpha = 0.6
    ) +
    ggplot2::geom_ribbon(
      ggplot2::aes_string(ymin = "rejected_died_q975", ymax = "rejected_died_q99"),
      fill = "yellow1",
      alpha = 0.6
    ) +
    ggplot2::geom_line(ggplot2::aes_string(y = "rejected_died_median")) +
    ggplot2::labs(title = "Outputs: Capacity-dependent deaths (per day)") +
    ggplot2::theme(
      axis.title.x = ggplot2::element_blank(),
      axis.title.y = ggplot2::element_blank(),
      plot.title = ggplot2::element_text(size = 11),
      axis.text.x = ggplot2::element_text(angle = 45, hjust = 1)
    ) +
    ggplot2::scale_x_date(date_breaks = "months", date_labels = "%b-%y")

  if(display){print(plot)}
  invisible(plot)
}


#' Plot cumulative total of admitted patients over simulation period
#'
#' @inheritParams plot_daily_hospitalisations
#' @inheritParams plot_cum_total_nonadmitted_surv
#' @param data_cum \code{data.frame}, cummalative output data.
#'
#' @return Plot of cumulative total addmitted patients over time. Invisibly returns the plot object.
#' @export
#'
#' @examples
plot_cum_admissions <- function(data_cum, display=TRUE) {
  #cumulative total of admitted patients over simulation period
  plot <- data_cum %>%
    ggplot2::ggplot(ggplot2::aes_string(x = "dates")) +
    #admitted
    ggplot2::geom_ribbon(ggplot2::aes_string(ymin = "admitted_q01", ymax = "admitted_q025"),
                         fill = "darkorange",
                         alpha = 0.6) +
    ggplot2::geom_ribbon(ggplot2::aes_string(ymin = "admitted_q025", ymax = "admitted_q05"),
                         fill = "darkorange1",
                         alpha = 0.6) +
    ggplot2::geom_ribbon(ggplot2::aes_string(ymin = "admitted_q05", ymax = "admitted_q15"),
                         fill = "darkorange2",
                         alpha = 0.6) +
    ggplot2::geom_ribbon(ggplot2::aes_string(ymin = "admitted_q15", ymax = "admitted_q30"),
                         fill = "darkorange3",
                         alpha = 0.6) +
    ggplot2::geom_ribbon(ggplot2::aes_string(ymin = "admitted_q30", ymax = "admitted_q70"),
                         fill = "darkorange4",
                         alpha = 0.6) +
    ggplot2::geom_ribbon(ggplot2::aes_string(ymin = "admitted_q70", ymax = "admitted_q85"),
                         fill = "darkorange3",
                         alpha = 0.6) +
    ggplot2::geom_ribbon(ggplot2::aes_string(ymin = "admitted_q85", ymax = "admitted_q95"),
                         fill = "darkorange2",
                         alpha = 0.6) +
    ggplot2::geom_ribbon(ggplot2::aes_string(ymin = "admitted_q95", ymax = "admitted_q975"),
                         fill = "darkorange1",
                         alpha = 0.6) +
    ggplot2::geom_ribbon(ggplot2::aes_string(ymin = "admitted_q975", ymax = "admitted_q99"),
                         fill = "darkorange",
                         alpha = 0.6) +
    ggplot2::geom_line(ggplot2::aes_string(y = "admitted_median")) +
    ggplot2::labs(title = "Outputs: Cumulative total for admitted patients") +
    ggplot2::scale_x_date(date_breaks = "months", date_labels = "%b-%y") +
    ggplot2::ylim(
      0,
      max(
        data_cum$admitted_q99,
        data_cum$rejected_died_q99,
        data_cum$rejected_survived_q99
      )
    ) +
    ggplot2::theme(
      axis.title.x = ggplot2::element_blank(),
      axis.title.y = ggplot2::element_blank(),
      plot.title = ggplot2::element_text(size = 11),
      axis.text.x = ggplot2::element_text(angle = 45, hjust = 1)
    )

  if(display){print(plot)}
  invisible(plot)
}


#' Plot cumulative total of patients who could not be admitted because of capacity constraints but did
#'   surivive, over simulation period
#'
#' @inheritParams plot_daily_hospitalisations
#' @inheritParams plot_cum_admissions
#'
#' @return Plot of those not admitted but survived over time. Invisibly returns the plot object.
#' @export
#'
#' @examples
plot_cum_total_nonadmitted_surv <- function(data_cum, display=TRUE) {
  #cumulative total of patients who could not be admitted because of capacity constraints but did
  #surivive, over simulation period
  plot <- data_cum %>%
    ggplot2::ggplot(ggplot2::aes_string(x = "dates")) +
    #rejected-survived
    ggplot2::geom_ribbon(
      ggplot2::aes_string(ymin = "rejected_survived_q01", ymax = "rejected_survived_q025"),
      fill = "chartreuse",
      alpha = 0.6
    ) +
    ggplot2::geom_ribbon(
      ggplot2::aes_string(ymin = "rejected_survived_q025", ymax = "rejected_survived_q05"),
      fill = "chartreuse1",
      alpha = 0.6
    ) +
    ggplot2::geom_ribbon(
      ggplot2::aes_string(ymin = "rejected_survived_q05", ymax = "rejected_survived_q15"),
      fill = "chartreuse2",
      alpha = 0.6
    ) +
    ggplot2::geom_ribbon(
      ggplot2::aes_string(ymin = "rejected_survived_q15", ymax = "rejected_survived_q30"),
      fill = "chartreuse3",
      alpha = 0.6
    ) +
    ggplot2::geom_ribbon(
      ggplot2::aes_string(ymin = "rejected_survived_q30", ymax = "rejected_survived_q70"),
      fill = "chartreuse4",
      alpha = 0.6
    ) +
    ggplot2::geom_ribbon(
      ggplot2::aes_string(ymin = "rejected_survived_q70", ymax = "rejected_survived_q85"),
      fill = "chartreuse3",
      alpha = 0.6
    ) +
    ggplot2::geom_ribbon(
      ggplot2::aes_string(ymin = "rejected_survived_q85", ymax = "rejected_survived_q95"),
      fill = "chartreuse2",
      alpha = 0.6
    ) +
    ggplot2::geom_ribbon(
      ggplot2::aes_string(ymin = "rejected_survived_q95", ymax = "rejected_survived_q975"),
      fill = "chartreuse1",
      alpha = 0.6
    ) +
    ggplot2::geom_ribbon(
      ggplot2::aes_string(ymin = "rejected_survived_q975", ymax = "rejected_survived_q99"),
      fill = "chartreuse",
      alpha = 0.6
    ) +
    ggplot2::geom_line(ggplot2::aes_string(y = "rejected_survived_median")) +
    ggplot2::labs(title = "Outputs: Cumulative total for rejected and survived") +
    ggplot2::scale_x_date(date_breaks = "months", date_labels = "%b-%y") +
    ggplot2::ylim(
      0,
      max(
        data_cum$admitted_q99,
        data_cum$rejected_died_q99,
        data_cum$rejected_survived_q99
      )
    ) +
    ggplot2::theme(
      axis.title.x = ggplot2::element_blank(),
      axis.title.y = ggplot2::element_blank(),
      plot.title = ggplot2::element_text(size = 11),
      axis.text.x = ggplot2::element_text(angle = 45, hjust = 1)
    )

  if(display){print(plot)}
  invisible(plot)
}

#' Plot cumulative total of patients who died as a result of insufficent capacity over simulation period
#'
#' @inheritParams plot_daily_hospitalisations
#' @inheritParams plot_cum_admissions
#'
#' @return Plot of cumulative patients dead due to lack of capacity. Invisibly returns the plot object.
#' @export
#'
#' @examples
plot_cum_died_inscap <- function(data_cum, display=TRUE) {
  #cumulative total of patients who died as a result of insufficent capacity over simulation period
  plot <- data_cum %>%
    ggplot2::ggplot(ggplot2::aes_string(x = "dates")) +
    #rejected-died
    ggplot2::geom_ribbon(
      ggplot2::aes_string(ymin = "rejected_died_q01", ymax = "rejected_died_q025"),
      fill = "yellow1",
      alpha = 0.6
    ) +
    ggplot2::geom_ribbon(
      ggplot2::aes_string(ymin = "rejected_died_q025", ymax = "rejected_died_q05"),
      fill = "orange2",
      alpha = 0.6
    ) +
    ggplot2::geom_ribbon(
      ggplot2::aes_string(ymin = "rejected_died_q05", ymax = "rejected_died_q15"),
      fill = "orangered3",
      alpha = 0.6
    ) +
    ggplot2::geom_ribbon(
      ggplot2::aes_string(ymin = "rejected_died_q15", ymax = "rejected_died_q30"),
      fill = "red3",
      alpha = 0.6
    ) +
    ggplot2::geom_ribbon(
      ggplot2::aes_string(ymin = "rejected_died_q30", ymax = "rejected_died_q70"),
      fill = "red2",
      alpha = 0.6
    ) +
    ggplot2::geom_ribbon(
      ggplot2::aes_string(ymin = "rejected_died_q70", ymax = "rejected_died_q85"),
      fill = "red3",
      alpha = 0.6
    ) +
    ggplot2::geom_ribbon(
      ggplot2::aes_string(ymin = "rejected_died_q85", ymax = "rejected_died_q95"),
      fill = "orangered3",
      alpha = 0.6
    ) +
    ggplot2::geom_ribbon(
      ggplot2::aes_string(ymin = "rejected_died_q95", ymax = "rejected_died_q975"),
      fill = "orange2",
      alpha = 0.6
    ) +
    ggplot2::geom_ribbon(
      ggplot2::aes_string(ymin = "rejected_died_q975", ymax = "rejected_died_q99"),
      fill = "yellow1",
      alpha = 0.6
    ) +
    ggplot2::geom_line(ggplot2::aes_string(y = "rejected_died_median")) +
    ggplot2::labs(title = "Outputs: Cumulative total for rejected and died") +
    ggplot2::scale_x_date(date_breaks = "months", date_labels = "%b-%y") +
    ggplot2::ylim(
      0,
      max(
        data_cum$admitted_q99,
        data_cum$rejected_died_q99,
        data_cum$rejected_survived_q99
      )
    ) +
    ggplot2::theme(
      axis.title.x = ggplot2::element_blank(),
      axis.title.y = ggplot2::element_blank(),
      plot.title = ggplot2::element_text(size = 11),
      axis.text.x = ggplot2::element_text(angle = 45, hjust = 1)
    )

  if(display){print(plot)}
  invisible(plot)
}

#' Hospitalisation time-series plotting function
#'
#' This is a convenience function to plot the hospitalisations over time for a given input dataset for the simulation
#' @inheritParams plot_daily_hospitalisations
#' @param case_df \code{data.frame}, cases admitted with two columns, \code{dates}; \code{str} dates, \code{hospitalisations}; \code{int} number of cases.
#' @format A data frame with 2 variables:
#' \describe{
#'   \item{dates}{\code{string}, date}
#'   \item{hospitalisations}{\code{int}, hopitalisation cases}
#' }
#'
#' @return Time-series plot of hospitalisations. Invisibly returns the plot object.
#' @export
#'
#' @examples
plot_case_scenario <- function(case_df, display = TRUE) {
  dates <- NULL # preallocate

  if(any(!c("hospitalisations", "dates") %in% names(case_df))) stop("Data does not contain the variables 'hospitalisations' and 'dates'")
  case_df <- case_df %>%
    dplyr::mutate(case_df, dates = lubridate::dmy(dates))
  p <- ggplot2::ggplot(case_df,
                       ggplot2::aes_string(x = "dates",
                                           y = "hospitalisations"))
  p <- p + ggplot2::geom_line()
  p <- p + ggplot2::labs(title = "Hospitlisations over time",
                         x = "Date")
  print(p)
}