fix plotting list instead of single sim

NAMESPACE

 # Generated by roxygen2: do not edit by hand
 
 S3method(as.data.frame,seiqhrf)
+S3method(plot,list)
 S3method(plot,seiqhrf)
-S3method(plot,sirplus)
 S3method(print,control.seiqhrf)
 S3method(print,init.seiqhrf)
 S3method(print,param.seiqhrf)
@@ -22,6 +22,7 @@ export(init_seiqhrf)
 export(init_status.icm)
 export(initialize.FUN)
 export(param_seiqhrf)
+export(plot_sirplus)
 export(plot_times)
 export(recovery.FUN)
 export(seiqhrf)

R/FIN_plot.R

 #' Plot simulation result 
 #'
-#' Function to extract timings, prevalence etc. from the simulation 
-#'
 #' @param x An seiqhrf object returned from function \code{\link{seiqhrf}}.
 #' @param method If "times", plot Duration frequency distributions.
 #'               If "weekly_local", plot local weekly estimates from simulation.
@@ -28,8 +26,6 @@
 #'        etc.
 #' @param ... Additional parameters
 #' 
-#' @return ggplot2 object
-#' 
 #' @importFrom tidyr pivot_longer
 #' @importFrom dplyr mutate
 #' @importFrom dplyr "%>%"
@@ -56,7 +52,7 @@ plot.seiqhrf <- function(x,
                          total_population = NULL, ...) {
     
     if(is.null(method)){
-        plot.sirplus(x, comp_remove = comp_remove,
+        plot_sirplus(x, comp_remove = comp_remove,
                      time_lim = time_lim,
                      ci = ci,
                      sep_compartments = sep_compartments,
@@ -86,6 +82,54 @@ plot.seiqhrf <- function(x,
     }
 }
 
+#' Plot simulation result given list
+#'
+#' @param x A list of seiqhrf objects returned from \code{\link{seiqhrf}}.
+#' @param return_df In effect only when method == "weekly", if TRUE returns 
+#'        also the dataframe used for plotting as well as the ggplot object.
+#' @param comp_remove Compartments to remove. Suggest c(s.num, r.num)
+#' @param time_lim Number of steps (days) to plot. 
+#' @param ci y/n to include 95% confidence intervals in sirplus plot.
+#' @param sep_compartments y/n use faceting to show each compartment in a 
+#'        separate plot, only works if plotting a single simulation.
+#' @param trans Y-axis transformation (e.g. log2, log10). Default = none. 
+#' @param known Dataframe with known compartment numbers to plot alongside
+#'        projections
+#' @param start_date Date for day 0. Default: ymd("2020-03-21"),
+#' @param x_axis Title for x-axis. Default: 'Days since beginning of epidemic'
+#' @param plot_title Title for whole plot. Default: 'SEIQHRF plot'
+#' @param ... Additional parameters
+#' 
+#' @importFrom tidyr pivot_longer
+#' @importFrom dplyr mutate
+#' @importFrom dplyr "%>%"
+#' @importFrom dplyr bind_rows
+#' @importFrom dplyr select
+#' @importFrom dplyr filter
+#' @import lubridate
+#' @import ggplot2
+#' @export
+plot.list <- function(x, comp_remove = "none",
+                      time_lim = 90,
+                      ci = 'y',
+                      sep_compartments = 'n',
+                      trans = 'na',
+                      known = NULL,
+                      start_date = ymd("2020-03-21"),
+                      x_axis = 'Days since beginning of epidemic',
+                      plot_title = 'SEIQHRF', ...){
+    
+    plot_sirplus(x, comp_remove = comp_remove,
+                 time_lim = time_lim,
+                 ci = ci,
+                 sep_compartments = sep_compartments,
+                 trans = trans,
+                 known = known,
+                 start_date = start_date,
+                 x_axis = x_axis,
+                 plot_title = plot_title)
+}
+
 
 #' Wrapper for primary sirplus plotting function
 #'
@@ -120,7 +164,7 @@ plot.seiqhrf <- function(x,
 #' @importFrom dplyr filter
 #' @import ggplot2
 #' @export
-plot.sirplus <- function(x,comp_remove = comp_remove,
+plot_sirplus <- function(x,comp_remove = comp_remove,
                          time_lim = time_lim,
                          ci = ci,
                          sep_compartments = sep_compartments,
@@ -312,7 +356,6 @@ get_ci <- function(x, plot_df){
     if(class(x) == "seiqhrf"){
         
         ci_info <- as.data.frame.list(summary.seiqhrf(x))
-        print(head(ci_info))
         ci_info <- ci_info %>% mutate(time = as.numeric(row.names(ci_info))) %>%
             pivot_longer(cols = -time, names_to = 'compartment',
                          values_to = 'mean') %>%

man/add_known.Rd

+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/FIN_plot.R
+\name{add_known}
+\alias{add_known}
+\title{Add known counts to sims dataframe for ggplot}
+\usage{
+add_known(plot_df, known = known, start_date = start_date)
+}
+\arguments{
+\item{known}{Dataframe with known compartment numbers to plot alongside
+projections}
+
+\item{x}{An seiqhrf object returned from function \code{\link{seiqhrf}}.}
+}
+\value{
+dataframe with known data added
+}
+\description{
+Add known counts to sims dataframe for ggplot
+}

man/format_sims.Rd

+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/FIN_plot.R
+\name{format_sims}
+\alias{format_sims}
+\title{Format seiqhrf objects into dataframe for ggplot}
+\usage{
+format_sims(x, time_lim = time_lim, start_date = start_date)
+}
+\arguments{
+\item{x}{An seiqhrf object returned from function \code{\link{seiqhrf}}.}
+
+\item{time_lim}{Number of steps (days) to plot.}
+
+\item{start_date}{Date for day 0. Default: ymd("2020-03-21"),}
+}
+\value{
+dataframe
+}
+\description{
+Format seiqhrf objects into dataframe for ggplot
+}

man/get_ci.Rd

+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/FIN_plot.R
+\name{get_ci}
+\alias{get_ci}
+\title{Get 95% confidence intervals}
+\usage{
+get_ci(x, plot_df)
+}
+\arguments{
+\item{x}{An seiqhrf object returned from function \code{\link{seiqhrf}}.}
+
+\item{known}{Dataframe with known compartment numbers to plot alongside
+projections}
+}
+\value{
+dataframe with CIs and sd added
+}
+\description{
+Get 95% confidence intervals
+}

man/plot.list.Rd

@@ -2,12 +2,14 @@
 % Please edit documentation in R/FIN_plot.R
 \name{plot.list}
 \alias{plot.list}
-\title{Plot simulation result}
+\title{Plot simulation result given list}
 \usage{
 \method{plot}{list}(
   x,
   comp_remove = "none",
-  time_lim = 100,
+  time_lim = 90,
+  ci = "y",
+  sep_compartments = "n",
   trans = "na",
   known = NULL,
   start_date = ymd("2020-03-21"),
@@ -17,12 +19,17 @@
 )
 }
 \arguments{
-\item{x}{A list of seiqhrf objects returned from function \code{\link{seiqhrf}}.}
+\item{x}{A list of seiqhrf objects returned from \code{\link{seiqhrf}}.}
 
 \item{comp_remove}{Compartments to remove. Suggest c(s.num, r.num)}
 
 \item{time_lim}{Number of steps (days) to plot.}
 
+\item{ci}{y/n to include 95% confidence intervals in sirplus plot.}
+
+\item{sep_compartments}{y/n use faceting to show each compartment in a 
+separate plot, only works if plotting a single simulation.}
+
 \item{trans}{Y-axis transformation (e.g. log2, log10). Default = none.}
 
 \item{known}{Dataframe with known compartment numbers to plot alongside
@@ -35,10 +42,10 @@ projections}
 \item{plot_title}{Title for whole plot. Default: 'SEIQHRF plot'}
 
 \item{...}{Additional parameters}
-}
-\value{
-ggplot2 object
+
+\item{return_df}{In effect only when method == "weekly", if TRUE returns 
+also the dataframe used for plotting as well as the ggplot object.}
 }
 \description{
-Function to extract timings, prevalence etc. from the simulation
+Plot simulation result given list
 }

man/plot.seiqhrf.Rd

@@ -8,7 +8,9 @@
   x,
   method = NULL,
   comp_remove = "none",
-  time_lim = 100,
+  time_lim = 90,
+  ci = "y",
+  sep_compartments = "n",
   trans = "na",
   known = NULL,
   start_date = ymd("2020-03-21"),
@@ -25,13 +27,18 @@
 \item{x}{An seiqhrf object returned from function \code{\link{seiqhrf}}.}
 
 \item{method}{If "times", plot Duration frequency distributions.
-If "weekly_local", plot local and weekly estimates from simulation.
-If NULL, plot individuals models or mutliple models for comparison.}
+If "weekly_local", plot local weekly estimates from simulation.
+If NULL, plot sirplus plots.}
 
 \item{comp_remove}{Compartments to remove. Suggest c(s.num, r.num)}
 
 \item{time_lim}{Number of steps (days) to plot.}
 
+\item{ci}{y/n to include 95% confidence intervals in sirplus plot.}
+
+\item{sep_compartments}{y/n use faceting to show each compartment in a 
+separate plot, only works if plotting a single simulation.}
+
 \item{trans}{Y-axis transformation (e.g. log2, log10). Default = none.}
 
 \item{known}{Dataframe with known compartment numbers to plot alongside
@@ -43,7 +50,8 @@ projections}
 
 \item{plot_title}{Title for whole plot. Default: 'SEIQHRF plot'}
 
-\item{return_df}{In effect only when method == "weekly", if TRUE returns also the dataframe used for plotting as well as the ggplot object}
+\item{return_df}{In effect only when method == "weekly", if TRUE returns 
+also the dataframe used for plotting as well as the ggplot object.}
 
 \item{market.share}{between 0 and 1, percentage of local hospital beds in 
 the simulated unit (e.g. state)}
@@ -57,9 +65,6 @@ etc.}
 
 \item{...}{Additional parameters}
 }
-\value{
-ggplot2 object
-}
 \description{
-Function to extract timings, prevalence etc. from the simulation
+Plot simulation result
 }

man/plot_sirplus.Rd

+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/FIN_plot.R
+\name{plot_sirplus}
+\alias{plot_sirplus}
+\title{Wrapper for primary sirplus plotting function}
+\usage{
+plot_sirplus(
+  x,
+  comp_remove = comp_remove,
+  time_lim = time_lim,
+  ci = ci,
+  sep_compartments = sep_compartments,
+  trans = trans,
+  known = known,
+  start_date = start_date,
+  x_axis = x_axis,
+  plot_title = plot_title,
+  ...
+)
+}
+\arguments{
+\item{x}{An seiqhrf object returned from function \code{\link{seiqhrf}}.}
+
+\item{comp_remove}{Compartments to remove. Suggest c(s.num, r.num)}
+
+\item{time_lim}{Number of steps (days) to plot.}
+
+\item{ci}{y/n to include 95% confidence intervals in sirplus plot.}
+
+\item{sep_compartments}{y/n use faceting to show each compartment in a 
+separate plot, only works if plotting a single simulation.}
+
+\item{trans}{Y-axis transformation (e.g. log2, log10). Default = none.}
+
+\item{known}{Dataframe with known compartment numbers to plot alongside
+projections}
+
+\item{start_date}{Date for day 0. Default: ymd("2020-03-21"),}
+
+\item{x_axis}{Title for x-axis. Default: 'Days since beginning of epidemic'}
+
+\item{plot_title}{Title for whole plot. Default: 'SEIQHRF plot'}
+
+\item{...}{Additional parameters}
+
+\item{method}{If "times", plot Duration frequency distributions.
+If NULL, plot sirplus plots.}
+}
+\value{
+ggplot2 object
+}
+\description{
+Flexible function to generate sirplus plots (i.e. compartment counts over 
+time). This function allows for plotting multiple experiments, viewing the 
+plots of different scales (e.g. log2), plotting compartments separately,
+adding 95% CIs, and plotting known data along side the simulations.
+}

man/plot_times.Rd

@@ -2,7 +2,7 @@
 % Please edit documentation in R/FIN_plot.R
 \name{plot_times}
 \alias{plot_times}
-\title{Plot simulation by time}
+\title{Plot compartment duration distributions}
 \usage{
 plot_times(sim)
 }
@@ -13,5 +13,6 @@ plot_times(sim)
 ggplot2 object
 }
 \description{
-Function to plot Duration frequency distributions.
+Function to plot Duration frequency distributions. If multiple simulations 
+were performed (nsim >1), durations from sims are appended to each other.
 }

vignettes/sirplus_intro_v2.Rmd

@@ -40,7 +40,8 @@ Again, for more information about these compartments and about the parameters
 that are associated with them see Tim Churches's [blog post](https://timchurches.github.io/blog/posts/2020-03-18-modelling-the-effects-of-public-health-interventions-on-covid-19-transmission-part-2/).
 
 ```{r, load package}
-library(sirplus)
+#library(sirplus)
+devtools::load_all(".")
 ```
 
 
@@ -101,8 +102,10 @@ plot(baseline_sim, "times")
 
 ### Plot baseline sirplus results
 
-To visualise your sirplus model, you can plot the change in prevalence (i.e. people) 
-over time in each compartment. 
+To visualise your sirplus model, you can plot the change in prevalence 
+(i.e. people) over time in each compartment. By default, this plotting 
+function shows the mean count across all simulations and the 95th quantile in 
+the ribbon. 
 
 ```{r viz prevalence}
 plot(baseline_sim, 
@@ -139,11 +142,10 @@ param <- param_seiqhrf(act.rate.e = act_rate_relax, act.rate.i = act_rate_relax
 sim_exp_relax <- seiqhrf(init, control, param)
 
 # Compare experiments 1 and 2 to the baseline simulation
-plot(list("Baseline" = baseline_sim, "Closures" = sim_exp, "Closures (2 mo)" = sim_exp_relax), 
-            start_date = lubridate::ymd("2020-01-01"),
-            comp_remove = c('s.num', 'r.num'),
-            plot_title = 'Closures Experiment')
-
+plot(list("Baseline" = baseline_sim, "Closures" = sim_exp, 
+          "Closures (2 mo)" = sim_exp_relax),
+     start_date = lubridate::ymd("2020-01-01"), 
+     comp_remove = c('s.num', 'r.num'), plot_title = 'Closures Experiment')
 ```
 
 From these results we see that this policy would likey reduce the peak number
@@ -151,7 +153,49 @@ of infections from 500 to 300 and would "flatten the curve" for infections and
 thus hospitalizations (the peak in number of infected people is lower, but the number of people infected declines more slowly after the peak of the epidemic than in the baseline model).
 
 
-## Visualize sirplus
+## Visualize sirplus: Advanced Plotting Options
+
+### A. Adding known data
+
+You can use the `plot()` function to include known compartment values alongside
+your simulations by supplying the known values as a dataframe to parameter 
+`known`. This can be helpful when wanting to know how well experiments are 
+simulating the epidemic progression to the point where you have data. 
+
+For example, if we know the hospitalization numbers have been growing 
+exponentially (by 0.2) at each step, we can see how that compares to our 
+experiments.
+
+```{r sirplus plot with known}
+known <- data.frame('time' = seq(1:30), 'h.num' =  seq(0, 10, by=0.2)[1:30])
+
+plot(list("Baseline" = baseline_sim, "Closures" = sim_exp, 
+          "Closures (2 mo)" = sim_exp_relax), time_lim = 45,
+     known = known,
+     start_date = lubridate::ymd("2020-01-01"), 
+     comp_remove = c('s.num', 'r.num', 'e.num', 'i.num', 'q.num'), 
+     plot_title = 'Closures Experiment')
+```
+From this, we can see that the closure for 2 months is slighly pessimistic, 
+while the closures simulation is close to the known hospitalization rate.
+
+
+### B. Plotting compartment separately
+
+You can use the `plot()` function to plot counts for each compartment 
+separately by specifying `sep_compartments` as 'y'. This can be useful when 
+visualizing compartments on different scales (e.g. i.num v f.num) without using 
+log transformations. 
+
+```{r sirplus plot by compartment, fig.height=8}
+plot(sim_exp, 
+     sep_compartments = 'y',
+     start_date = lubridate::ymd("2020-01-01"), 
+     plot_title = 'Closures Experiment')
+```
+
+### C. Weekly hospitalization numbers
+
 Use the `plot()` function with parameter `method = "weekly_local"` to extract and plot the number of weekly expected patients in a hospital. If `return_df` is `TRUE`, the function returns a dataframe that generates the plot, otherwise, it returns a ggplot object. 
 
 This function takes the following input: