fixed hosp plots

R/FIN_plot.R

@@ -30,7 +30,6 @@
 #' @importFrom dplyr mutate
 #' @importFrom dplyr "%>%"
 #' @importFrom dplyr bind_rows
-#' @importFrom dplyr select
 #' @importFrom dplyr filter
 #' @import lubridate
 #' @import ggplot2
@@ -72,7 +71,7 @@ plot.seiqhrf <- function(x, method = NULL,
             ret <- get_weekly_local(x, market.share = market.share,
                                     icu_percent = icu_percent, 
                                     start_date = start_date,
-                                    time_limit = time_lim,
+                                    time_lim = time_lim,
                                     total_population = total_population)
              if(return_df){
                  return(ret) 
@@ -102,7 +101,6 @@ plot.seiqhrf <- function(x, method = NULL,
 #' @importFrom dplyr mutate
 #' @importFrom dplyr "%>%"
 #' @importFrom dplyr bind_rows
-#' @importFrom dplyr select
 #' @importFrom dplyr filter
 #' @import lubridate
 #' @import ggplot2
@@ -157,7 +155,6 @@ plot.list <- function(x, comp_remove = "none",
 #' @importFrom dplyr mutate
 #' @importFrom dplyr "%>%"
 #' @importFrom dplyr bind_rows
-#' @importFrom dplyr select
 #' @importFrom dplyr filter
 #' @import ggplot2
 #' @export
@@ -293,6 +290,92 @@ plot_times <- function(sim) {
 }
 
 
+#' Extract and plot information of local and weekly estimates from simulation 
+#'
+#' @param sim An \code{seiqhrf} object returned by \link{simulate_seiqhrf}. 
+#' @param market.share between 0 and 1, percentage of local hospital beds in 
+#'        the simulated unit (e.g. state)
+#' @param icu_percent between 0 and 1, percentage of patients that should go to 
+#'        ICU among the ones that need hospitalization
+#' @param start_date Epidemic start date. Default is 'na', if not provided will 
+#'        plot week numbers, if provided will plot the first day (Sunday) of the
+#'        week.
+#' @param time_limit Number of days to include. Default = 90.
+#' @param total_population True population size, needed only if simulation size 
+#'        is smaller than the true population size due to computational cost 
+#'        etc.
+#' 
+#' @return 
+#' \itemize{
+#' \item \code{plot:} A \code{ggplot} object, bar charts of count of patients 
+#'              requiring hospitalization and ICU respectively
+#' \item \code{result:} A dataframe
+#'    \itemize{\item \code{week:}  week number from input \code{sim},
+#'             \item \code{hosp:} the number of patients that require hospitalization locally,
+#'             \item \code{icu:} the number of patients that require ICU locally. }
+#
+#' }
+#' 
+#' @importFrom tidyr pivot_wider
+#' @export
+get_weekly_local <- function(sim, 
+                             market.share = 0.04,
+                             icu_percent = .1, 
+                             start_date = ymd("2020-03-21"),
+                             time_lim = 90,
+                             total_population = NULL){
+
+    # Get h.num and 95% quantile CIs
+    sim_mean <- as.data.frame(sim, out = "mean")
+    ci_info <- as.data.frame.list(summary.seiqhrf(sim)$h.num)
+    hosp <- data.frame('h.num' = sim_mean$h.num, 'ci5' = ci_info$qntCI.1,
+                       'ci95' = ci_info$qntCI.2)
+    hosp[is.na(hosp)] <- 0
+    hosp <- hosp[1: time_lim, ]
+    hosp$date <- start_date + as.numeric(row.names(hosp))
+    
+    # Scale for population size and hospital market share if needed
+    if(!is.null(total_population)){
+        if(total_population < max(sim_mean$s.num)) 
+            stop("total Population should be larger than simulated size")
+        cat("Scalling w.r.t total population")
+        hosp <- hosp*total_population/max(sim_mean$s.num)
+    } 
+    
+    if(market.share < 0 || market.share > 1) stop("Market share has to be between 
+                                                0 and 1")
+    if(icu_percent < 0 || icu_percent > 1) stop("ICU percentage has to be between
+                                              0 and 1")
+    
+    # Get weekly sums & calculate projected icu numbers
+    hosp.wk <- hosp %>% group_by(yr_wk = floor_date(date, "1 week")) %>% 
+        summarise(h.num=sum(h.num), h.ci5=sum(ci5), h.ci95=sum(ci95)) %>%
+        mutate(icu.num = h.num * icu_percent,
+               icu.ci5 = h.ci5 * icu_percent,
+               icu.ci95 = h.ci95 * icu_percent) %>%
+        mutate(h.num = h.num - icu.num,
+               h.ci5 = h.ci5 - icu.ci5,
+               h.ci95 = h.ci95 - icu.ci95) 
+    
+    # Make long format for ggplot
+    hosp.wk2 <- hosp.wk %>% pivot_longer(-yr_wk, names_to = c('type', 'metric'), 
+                     values_to = 'val',
+                     names_pattern = '(h|icu).(num|ci5|ci95)') %>%
+        pivot_wider(names_from = metric, values_from = val)
+
+    p <- ggplot(hosp.wk2, aes(x = yr_wk, y = num, fill = type)) +
+        geom_bar(stat="identity") + theme_bw() +
+        scale_x_date(date_breaks = "1 week", date_labels = "%y-%m-%d") +
+        labs(y="Weekly Cumulative Count", x = "Week (Monday)") +
+        geom_errorbar(aes(ymin=ci5, ymax=ci95), size=0.5, width=.5) +
+        scale_fill_discrete(name = "Type", labels = c("General", "ICU"))
+    
+    return(list("plot" = p, "result" = hosp.wk))
+    
+}
+
+
+
 #' Format seiqhrf objects into dataframe for ggplot
 #'
 #' @param x An seiqhrf object returned from function \code{\link{seiqhrf}}.

man/get_weekly_local.Rd

 % Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/get_weekly_local.R
+% Please edit documentation in R/FIN_plot.R
 \name{get_weekly_local}
 \alias{get_weekly_local}
-\title{Extract information of local and weekly estimates from simulation}
+\title{Extract and plot information of local and weekly estimates from simulation}
 \usage{
 get_weekly_local(
   sim,
   market.share = 0.04,
   icu_percent = 0.1,
-  start_date = "na",
-  time_limit = 90,
+  start_date = ymd("2020-03-21"),
+  time_lim = 90,
   total_population = NULL
 )
 }
@@ -26,11 +26,11 @@ ICU among the ones that need hospitalization}
 plot week numbers, if provided will plot the first day (Sunday) of the
 week.}
 
-\item{time_limit}{Number of days to include. Default = 90.}
-
 \item{total_population}{True population size, needed only if simulation size 
 is smaller than the true population size due to computational cost 
 etc.}
+
+\item{time_limit}{Number of days to include. Default = 90.}
 }
 \value{
 \itemize{
@@ -43,5 +43,5 @@ etc.}
 }
 }
 \description{
-Extract information of local and weekly estimates from simulation
+Extract and plot information of local and weekly estimates from simulation
 }