From 0ec19e814abf4d45e9c7b52f38dde780eb43f614 Mon Sep 17 00:00:00 2001
From: pqiao29 <pqiao@student.unimelb.edu.au>
Date: Mon, 13 Apr 2020 16:33:46 +1000
Subject: [PATCH] S3 methods in progress: summary, plot (times, models)
---
NAMESPACE | 6 +-
R/{results-parse.R => FIN_plot.R} | 212 ++++++++++++------
R/{summary.R => FIN_summary.R} | 5 +-
R/{internal_update.R => INTN_update.R} | 0
R/PREP_init.R | 2 +-
R/original_plot_models.R | 123 ++++++++++
R/{RUN_simulation.R => original_simulation.R} | 5 +-
man/get_times.Rd | 17 --
man/init_seiqhrf.Rd | 4 +-
man/plot.list.Rd | 44 ++++
man/plot.seiqhrf.Rd | 48 ++++
man/plot_models.Rd | 6 +-
man/plot_times.Rd | 8 +-
man/simulate_seiqhrf.Rd | 7 +-
...{summary_seiqhrf.Rd => summary.seiqhrf.Rd} | 12 +-
tests/testthat/.DS_Store | Bin 6148 -> 6148 bytes
.../{test-seiqhrf.R => test-all-modules.R} | 2 +-
tests/testthat/test-plot-model.R | 49 ++++
tests/testthat/test-seiqhrf-simulate.R | 32 +++
vignettes/sirplus_intro.Rmd | 3 +-
vignettes/sirplus_intro_v2.Rmd | 82 ++++++-
21 files changed, 549 insertions(+), 118 deletions(-)
rename R/{results-parse.R => FIN_plot.R} (63%)
rename R/{summary.R => FIN_summary.R} (93%)
rename R/{internal_update.R => INTN_update.R} (100%)
create mode 100644 R/original_plot_models.R
rename R/{RUN_simulation.R => original_simulation.R} (99%)
delete mode 100644 man/get_times.Rd
create mode 100644 man/plot.list.Rd
create mode 100644 man/plot.seiqhrf.Rd
rename man/{summary_seiqhrf.Rd => summary.seiqhrf.Rd} (75%)
rename tests/testthat/{test-seiqhrf.R => test-all-modules.R} (91%)
create mode 100644 tests/testthat/test-plot-model.R
create mode 100644 tests/testthat/test-seiqhrf-simulate.R
diff --git a/NAMESPACE b/NAMESPACE
index 6494056..6aeed64 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -1,27 +1,27 @@
# Generated by roxygen2: do not edit by hand
S3method(as.data.frame,seiqhrf)
+S3method(plot,list)
+S3method(plot,seiqhrf)
S3method(print,control.seiqhrf)
S3method(print,init.seiqhrf)
S3method(print,param.seiqhrf)
S3method(print,seiqhrf)
+S3method(summary,seiqhrf)
export(arrivals.FUN)
export(control_seiqhrf)
export(cum_discr_si)
export(departures.FUN)
export(get_prev.FUN)
-export(get_times)
export(infection.FUN)
export(init_seiqhrf)
export(init_status.icm)
export(initialize.FUN)
export(param_seiqhrf)
-export(plot_models)
export(plot_times)
export(recovery.FUN)
export(seiqhrf)
export(simulate_seiqhrf)
-export(summary_seiqhrf)
export(vary_param)
import(dplyr)
import(ggplot2)
diff --git a/R/results-parse.R b/R/FIN_plot.R
similarity index 63%
rename from R/results-parse.R
rename to R/FIN_plot.R
index 83fcbfe..cd1d0c4 100644
--- a/R/results-parse.R
+++ b/R/FIN_plot.R
@@ -1,10 +1,75 @@
-#' Get data by time
+#' Plot simulation result
#'
-#' Function to extract timings and assemble results in a dataframe
+#' Function to extract timings, prevalence etc. from the simulation
#'
-#' @param simulate_results results from `simulte()` function.
+#' @param x An seiqhrf object returned from function \code{\link{seiqhrf}}.
+#' @param method If "times", plot Duration frequency distributions.
+#' If NULL, plot individuals models or mutliple models for comparison.
+#' @param comp_remove Compartments to remove. Suggest c(s.num, r.num)
+#' @param time_lim Number of steps (days) to plot.
+#' @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
+#'
+#' @return ggplot2 object
+#'
+#' @importFrom tidyr pivot_longer
+#' @importFrom dplyr mutate
+#' @importFrom dplyr "%>%"
+#' @importFrom dplyr bind_rows
+#' @importFrom dplyr select
+#' @importFrom dplyr filter
+#' @import ggplot2
+#' @export
+plot.seiqhrf <- function(x, method = NULL, comp_remove = "none",
+ time_lim = 100,
+ trans = 'na',
+ known = NULL,
+ start_date = ymd("2020-03-21"),
+ x_axis = 'Days since beginning of epidemic',
+ plot_title = 'SEIQHRF', ...) {
+
+ if(is.null(method)){
+
+ plot_model(x,
+ comp_remove = comp_remove,
+ time_lim = time_lim,
+ trans = trans,
+ known = known,
+ start_date = start_date,
+ x_axis = x_axis,
+ plot_title = plot_title)
+ }else{
+ if(method == "times"){
+ if(!inherits(x, "seiqhrf")) stop("If method == times, x needs to be an seiqhrf object")
+ plot_times(x)
+ }
+ }
+
+}
+
+
+
+#' Plot simulation result
+#'
+#' Function to extract timings, prevalence etc. from the simulation
+#'
+#' @param x A list of seiqhrf objects returned from function \code{\link{seiqhrf}}.
+#' @param comp_remove Compartments to remove. Suggest c(s.num, r.num)
+#' @param time_lim Number of steps (days) to plot.
+#' @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
#'
-#' @return dataframe with sim.n, period type, and duration
+#' @return ggplot2 object
#'
#' @importFrom tidyr pivot_longer
#' @importFrom dplyr mutate
@@ -12,10 +77,44 @@
#' @importFrom dplyr bind_rows
#' @importFrom dplyr select
#' @importFrom dplyr filter
+#' @import ggplot2
#' @export
-get_times <- function(simulate_results) {
+plot.list <- function(x, comp_remove = "none",
+ time_lim = 100,
+ trans = 'na',
+ known = NULL,
+ start_date = ymd("2020-03-21"),
+ x_axis = 'Days since beginning of epidemic',
+ plot_title = 'SEIQHRF', ...){
- sim <- simulate_results$sim
+ plot_model(x,
+ comp_remove = comp_remove,
+ time_lim = time_lim,
+ trans = trans,
+ known = known,
+ start_date = start_date,
+ x_axis = x_axis,
+ plot_title = plot_title)
+}
+
+#' Plot simulation by time
+#'
+#' Function to plot Duration frequency distributions.
+#'
+#' @param sim An seiqhrf object returned from function \code{\link{seiqhrf}}.
+#'
+#' @return ggplot2 object
+#'
+#' @import ggplot2
+#' @importFrom tidyr pivot_longer
+#' @importFrom dplyr mutate
+#' @importFrom dplyr "%>%"
+#' @importFrom dplyr bind_rows
+#' @importFrom dplyr select
+#' @importFrom dplyr filter
+#'
+#' @export
+plot_times <- function(sim) {
for (s in 1:sim$control$nsims) {
if (s == 1) {
@@ -52,80 +151,40 @@ get_times <- function(simulate_results) {
"Hosp care required duration",
"Survival time case fatalities"),
ordered = TRUE))
-
- return(times)
-}
-
-
-#' Plot times
-#'
-#' Function to plot Duration frequency distributions.
-#'
-#' @param times results from `get_times()` function.
-#'
-#' @return ggplot2 object
-#' @import ggplot2
-#' @export
-plot_times <- function(times){
+
times %>% filter(duration <= 30) %>% ggplot(aes(x = duration)) +
geom_bar() + facet_grid(period_type ~ ., scales = "free_y") +
labs(title = "Duration frequency distributions",
subtitle = "Baseline simulation")
}
-#' Plot models
-#'
-#' Function to plot individuals models or mutliple models for comparison.
-#'
-#' @param sims Single or list of sims to plot
-#' @param sim_id String or list of strings to name each facet
-#' @param comp_remove Compartments to remove. Suggest c(s.num, r.num)
-#' @param time_lim Number of steps (days) to plot.
-#' @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: 'ICM plot'
-#'
-#' @return ggplot2 object
-#'
-#' @import dplyr
-#' @importFrom tidyr pivot_longer
-#'
-#' @export
-plot_models <- function(sims = baseline_sim,
- sim_id = "baseline",
- comp_remove = "none",
- time_lim = 100,
- trans = 'na',
- known = NULL,
- start_date = ymd("2020-03-21"),
- x_axis = 'Days since beginning of epidemic',
- plot_title = 'ICM plot'){
-
- # Define a standard set of colours to represent compartments
- comps <- c("s.num", "e.num", "i.num", "q.num", "h.num", "r.num", "f.num")
- compcols <- c(s.num = "#4477AA", e.num = "#66CCEE", i.num = "#CCBB44",
- q.num = "#AA3377", h.num = "#EE6677", r.num = "#228833",
- f.num = "#BBBBBB")
- complabels <- c(s.num = "S: Susceptible", e.num = "E: Asymptomatic",
- i.num = "I: Infected", q.num = "Q: Self-isolated",
- h.num = "H: Hospitalized", r.num = "R: Recovered",
- f.num = "F: Case Fatalities")
+
+
+# Internal
+plot_model <- function(sims, comp_remove, time_lim, trans, known, start_date, x_axis, plot_title){
# Merge models to plot together
- for (i in (1: length(sim_id))) {
- if (i == 1) {
- plot_df <- sims[i*2]$df
- plot_df <- plot_df %>% mutate(experiment = sim_id[i])
- } else{
- tmp_df <- sims[i*2]$df
- tmp_df <- tmp_df %>% mutate(experiment = sim_id[i])
- plot_df <- plot_df %>% bind_rows(plot_df, tmp_df)
+ if(class(sims) == "seiqhrf"){
+ sim_id <- "Baseline"
+ plot_df <- as.data.frame(sims, out = "mean")
+ plot_df <- plot_df %>% mutate(experiment = sim_id)
+ }else{
+
+ sim_id <- names(sims)
+ if(is.null(sim_id)) stop("Please assign a name to each element in sims")
+
+ plot_df <- as.data.frame(sims[[1]], out = "mean")
+ plot_df <- plot_df %>% mutate(experiment = sim_id[1])
+ if(length(sim_id) > 1){
+ for (i in (2:length(sim_id))) {
+ tmp_df <- as.data.frame(sims[[i]], out = "mean")
+ tmp_df <- tmp_df %>% mutate(experiment = sim_id[i])
+ plot_df <- plot_df %>% bind_rows(plot_df, tmp_df)
+ }
}
}
+
plot_df <- plot_df %>% filter(time <= time_lim) %>%
pivot_longer(-c(time, experiment),
names_to = "compartment",
@@ -134,6 +193,17 @@ plot_models <- function(sims = baseline_sim,
plot_df$Date <- start_date + plot_df$time
plot_df$time <- NULL
+
+ # Define a standard set of colours to represent compartments
+ comps <- c("s.num", "e.num", "i.num", "q.num", "h.num", "r.num", "f.num")
+ compcols <- c(s.num = "#4477AA", e.num = "#66CCEE", i.num = "#CCBB44",
+ q.num = "#AA3377", h.num = "#EE6677", r.num = "#228833",
+ f.num = "#BBBBBB")
+ complabels <- c(s.num = "S: Susceptible", e.num = "E: Asymptomatic",
+ i.num = "I: Infected", q.num = "Q: Self-isolated",
+ h.num = "H: Hospitalized", r.num = "R: Recovered",
+ f.num = "F: Case Fatalities")
+
# Add known data
if(is.data.frame(known)){
exps <- unique(plot_df$experiment)
@@ -146,7 +216,7 @@ plot_models <- function(sims = baseline_sim,
# Filter compartments
comp_plot <- setdiff(comps, comp_remove)
plot_df <- plot_df %>% filter(compartment %in% c(comp_plot))
-
+
reo_exp <- function(x) { factor(x, levels = sim_id)}
@@ -197,3 +267,5 @@ plot_models <- function(sims = baseline_sim,
)
}
}
+
+
diff --git a/R/summary.R b/R/FIN_summary.R
similarity index 93%
rename from R/summary.R
rename to R/FIN_summary.R
index dfc0b2b..b8902b9 100644
--- a/R/summary.R
+++ b/R/FIN_summary.R
@@ -2,7 +2,8 @@
#'
#' Function to extract mean and sdandard deviation of all compartments for SEIQHRF model
#'
-#' @param object ICM parameters.
+#' @param object seiqhrf object, returned by \code{\link{seiqhrf}}.
+#' @param ... Additional parameters.
#'
#' @return A list of compartments, each compartment is a list of three:
#' \itemize{
@@ -14,7 +15,7 @@
#' where t is the total number of time points in the simulation (i.e. object$nstep).
#'
#' @export
-summary_seiqhrf <- function(object){
+summary.seiqhrf <- function(object, ...){
nsteps <- object$control$nsteps
prev_names <- c("s.num", "e.num", "i.num", "q.num", "h.num", "r.num", "f.num")
diff --git a/R/internal_update.R b/R/INTN_update.R
similarity index 100%
rename from R/internal_update.R
rename to R/INTN_update.R
diff --git a/R/PREP_init.R b/R/PREP_init.R
index 7f8147b..4c311b0 100644
--- a/R/PREP_init.R
+++ b/R/PREP_init.R
@@ -37,7 +37,7 @@
#' @export
#'
-init_seiqhrf <- function(s.num, e.num=0, i.num = 0, q.num=0, h.num=0, r.num = 0, f.num = 0, ...) {
+init_seiqhrf <- function(s.num = 9997, e.num=0, i.num = 3, q.num=0, h.num=0, r.num = 0, f.num = 0, ...) {
p <- list()
formal.args <- formals(sys.function())
diff --git a/R/original_plot_models.R b/R/original_plot_models.R
new file mode 100644
index 0000000..6842218
--- /dev/null
+++ b/R/original_plot_models.R
@@ -0,0 +1,123 @@
+#' Plot models
+#'
+#' Function to plot individuals models or mutliple models for comparison.
+#'
+#' @param sims Single or list of sims to plot
+#' @param sim_id String or list of strings to name each facet
+#' @param comp_remove Compartments to remove. Suggest c(s.num, r.num)
+#' @param time_lim Number of steps (days) to plot.
+#' @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'
+#'
+#' @return ggplot2 object
+#'
+#' @import dplyr
+#' @importFrom tidyr pivot_longer
+
+plot_models <- function(sims = baseline_sim,
+ sim_id = "baseline",
+ comp_remove = "none",
+ time_lim = 100,
+ trans = 'na',
+ known = NULL,
+ start_date = ymd("2020-03-21"),
+ x_axis = 'Days since beginning of epidemic',
+ plot_title = 'SEIQHRF'){
+
+ # Define a standard set of colours to represent compartments
+ comps <- c("s.num", "e.num", "i.num", "q.num", "h.num", "r.num", "f.num")
+ compcols <- c(s.num = "#4477AA", e.num = "#66CCEE", i.num = "#CCBB44",
+ q.num = "#AA3377", h.num = "#EE6677", r.num = "#228833",
+ f.num = "#BBBBBB")
+ complabels <- c(s.num = "S: Susceptible", e.num = "E: Asymptomatic",
+ i.num = "I: Infected", q.num = "Q: Self-isolated",
+ h.num = "H: Hospitalized", r.num = "R: Recovered",
+ f.num = "F: Case Fatalities")
+
+ # Merge models to plot together
+ for (i in (1: length(sim_id))) {
+ if (i == 1) {
+ plot_df <- sims[i*2]$df
+ plot_df <- plot_df %>% mutate(experiment = sim_id[i])
+ } else{
+ tmp_df <- sims[i*2]$df
+ tmp_df <- tmp_df %>% mutate(experiment = sim_id[i])
+ plot_df <- plot_df %>% bind_rows(plot_df, tmp_df)
+ }
+ }
+
+ plot_df <- plot_df %>% filter(time <= time_lim) %>%
+ pivot_longer(-c(time, experiment),
+ names_to = "compartment",
+ values_to = "count")
+ plot_df$sim <- 'sim'
+ plot_df$Date <- start_date + plot_df$time
+ plot_df$time <- NULL
+
+ # Add known data
+ if(is.data.frame(known)){
+ exps <- unique(plot_df$experiment)
+ for(i in exps){
+ known$experiment <- i
+ plot_df <- rbind(plot_df, known)
+ }
+ }
+
+ # Filter compartments
+ comp_plot <- setdiff(comps, comp_remove)
+ plot_df <- plot_df %>% filter(compartment %in% c(comp_plot))
+
+ reo_exp <- function(x) { factor(x, levels = sim_id)}
+
+
+ # Plot single model
+ if(trans == "na"){
+ if(length(sim_id) == 1){
+ plot_df %>% ggplot(aes(x = Date, y = count+1, colour = compartment,
+ linetype = sim)) +
+ geom_line(size = 1.5, alpha = 0.8) +
+ scale_x_date(date_breaks = "1 week", date_labels = "%m-%d") +
+ scale_colour_manual(values = compcols, labels = complabels) +
+ labs(title = plot_title, x = x_axis, y = "Prevalence") +
+ theme_bw() +
+ theme(axis.text.x = element_text(angle = 90))
+ } else (
+ plot_df %>% ggplot(aes(x = Date, y = count, colour = compartment,
+ linetype = sim)) +
+ facet_grid(reo_exp(experiment) ~ ., scale = 'free') +
+ scale_x_date(date_breaks = "1 week", date_labels = "%m-%d") +
+ geom_line(size = 1.5, alpha = 0.8) +
+ scale_colour_manual(values = compcols, labels = complabels) +
+ labs(title = plot_title, x = x_axis, y = "Prevalence") +
+ theme_bw() +
+ theme(axis.text.x = element_text(angle = 90))
+ )
+ }else{
+ if(length(sim_id) == 1){
+ plot_df %>% ggplot(aes(x = Date, y = count+1, colour = compartment,
+ linetype = sim)) +
+ scale_y_continuous(trans = trans) +
+ scale_x_date(date_breaks = "1 week", date_labels = "%m-%d") +
+ geom_line(size = 1.5, alpha = 0.8) +
+ scale_colour_manual(values = compcols, labels = complabels) +
+ labs(title = plot_title, x = x_axis, y = "Prevalence") +
+ theme_bw() +
+ theme(axis.text.x = element_text(angle = 90))
+ } else (
+ plot_df %>% ggplot(aes(x = Date, y = count+1, colour = compartment,
+ linetype = sim)) +
+ facet_grid(reo_exp(experiment) ~ ., scale = 'free') +
+ scale_y_continuous(trans = trans) +
+ scale_x_date(date_breaks = "1 week", date_labels = "%m-%d") +
+ geom_line(size = 1.5, alpha = 0.8) +
+ scale_colour_manual(values = compcols, labels = complabels) +
+ labs(title = plot_title, x = x_axis, y = "Prevalence") +
+ theme_bw() +
+ theme(axis.text.x = element_text(angle = 90))
+ )
+ }
+}
diff --git a/R/RUN_simulation.R b/R/original_simulation.R
similarity index 99%
rename from R/RUN_simulation.R
rename to R/original_simulation.R
index 07274ae..1e6783a 100644
--- a/R/RUN_simulation.R
+++ b/R/original_simulation.R
@@ -209,6 +209,7 @@
#' to model emigration as well as deaths.
#' @param out Summary function for the simulation runs. median is
#' also available, or percentiles, see the EpiModel documentation.
+#' @param seed Random see for checking consistency.
#'
#' @return list with simulation and simulation dataframe
#' @export
@@ -276,7 +277,7 @@ simulate_seiqhrf <- function(type = "SEIQHRF",
dq.rate = (7/365)/1000,
dh.rate = (20/365)/1000,
dr.rate = (7/365)/1000,
- out="mean"
+ out="mean", seed = NULL
) {
control <- control_seiqhrf(type = type,
@@ -343,7 +344,7 @@ simulate_seiqhrf <- function(type = "SEIQHRF",
dh.rate = dh.rate,
dr.rate = dr.rate)
- sim <- seiqhrf(param = param, init = init, control = control)
+ sim <- seiqhrf(param = param, init = init, control = control, seed = seed)
sim_df <- as.data.frame(sim, out=out)
class(sim) <- "icm"
diff --git a/man/get_times.Rd b/man/get_times.Rd
deleted file mode 100644
index e1b6676..0000000
--- a/man/get_times.Rd
+++ /dev/null
@@ -1,17 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/results-parse.R
-\name{get_times}
-\alias{get_times}
-\title{Get data by time}
-\usage{
-get_times(simulate_results)
-}
-\arguments{
-\item{simulate_results}{results from `simulte()` function.}
-}
-\value{
-dataframe with sim.n, period type, and duration
-}
-\description{
-Function to extract timings and assemble results in a dataframe
-}
diff --git a/man/init_seiqhrf.Rd b/man/init_seiqhrf.Rd
index 435306a..1676bed 100644
--- a/man/init_seiqhrf.Rd
+++ b/man/init_seiqhrf.Rd
@@ -5,9 +5,9 @@
\title{Initial Conditions for Stochastic Individual Contact Models}
\usage{
init_seiqhrf(
- s.num,
+ s.num = 9997,
e.num = 0,
- i.num = 0,
+ i.num = 3,
q.num = 0,
h.num = 0,
r.num = 0,
diff --git a/man/plot.list.Rd b/man/plot.list.Rd
new file mode 100644
index 0000000..820a51f
--- /dev/null
+++ b/man/plot.list.Rd
@@ -0,0 +1,44 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/FIN_plot.R
+\name{plot.list}
+\alias{plot.list}
+\title{Plot simulation result}
+\usage{
+\method{plot}{list}(
+ x,
+ comp_remove = "none",
+ time_lim = 100,
+ trans = "na",
+ known = NULL,
+ start_date = ymd("2020-03-21"),
+ x_axis = "Days since beginning of epidemic",
+ plot_title = "SEIQHRF",
+ ...
+)
+}
+\arguments{
+\item{x}{A list of seiqhrf objects 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{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}
+}
+\value{
+ggplot2 object
+}
+\description{
+Function to extract timings, prevalence etc. from the simulation
+}
diff --git a/man/plot.seiqhrf.Rd b/man/plot.seiqhrf.Rd
new file mode 100644
index 0000000..86009a2
--- /dev/null
+++ b/man/plot.seiqhrf.Rd
@@ -0,0 +1,48 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/FIN_plot.R
+\name{plot.seiqhrf}
+\alias{plot.seiqhrf}
+\title{Plot simulation result}
+\usage{
+\method{plot}{seiqhrf}(
+ x,
+ method = NULL,
+ comp_remove = "none",
+ time_lim = 100,
+ trans = "na",
+ known = NULL,
+ start_date = ymd("2020-03-21"),
+ x_axis = "Days since beginning of epidemic",
+ plot_title = "SEIQHRF",
+ ...
+)
+}
+\arguments{
+\item{x}{An seiqhrf object returned from function \code{\link{seiqhrf}}.}
+
+\item{method}{If "times", plot Duration frequency distributions.
+If NULL, plot individuals models or mutliple models for comparison.}
+
+\item{comp_remove}{Compartments to remove. Suggest c(s.num, r.num)}
+
+\item{time_lim}{Number of steps (days) to plot.}
+
+\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}
+}
+\value{
+ggplot2 object
+}
+\description{
+Function to extract timings, prevalence etc. from the simulation
+}
diff --git a/man/plot_models.Rd b/man/plot_models.Rd
index a6a19de..c158b0e 100644
--- a/man/plot_models.Rd
+++ b/man/plot_models.Rd
@@ -1,5 +1,5 @@
% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/results-parse.R
+% Please edit documentation in R/original_plot_models.R
\name{plot_models}
\alias{plot_models}
\title{Plot models}
@@ -13,7 +13,7 @@ plot_models(
known = NULL,
start_date = ymd("2020-03-21"),
x_axis = "Days since beginning of epidemic",
- plot_title = "ICM plot"
+ plot_title = "SEIQHRF"
)
}
\arguments{
@@ -34,7 +34,7 @@ projections}
\item{x_axis}{Title for x-axis. Default: 'Days since beginning of epidemic'}
-\item{plot_title}{Title for whole plot. Default: 'ICM plot'}
+\item{plot_title}{Title for whole plot. Default: 'SEIQHRF plot'}
}
\value{
ggplot2 object
diff --git a/man/plot_times.Rd b/man/plot_times.Rd
index cc2a1de..ca596ed 100644
--- a/man/plot_times.Rd
+++ b/man/plot_times.Rd
@@ -1,13 +1,13 @@
% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/results-parse.R
+% Please edit documentation in R/FIN_plot.R
\name{plot_times}
\alias{plot_times}
-\title{Plot times}
+\title{Plot simulation by time}
\usage{
-plot_times(times)
+plot_times(sim)
}
\arguments{
-\item{times}{results from `get_times()` function.}
+\item{sim}{An seiqhrf object returned from function \code{\link{seiqhrf}}.}
}
\value{
ggplot2 object
diff --git a/man/simulate_seiqhrf.Rd b/man/simulate_seiqhrf.Rd
index 4e4ae0e..782e580 100644
--- a/man/simulate_seiqhrf.Rd
+++ b/man/simulate_seiqhrf.Rd
@@ -1,5 +1,5 @@
% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/RUN_simulation.R
+% Please edit documentation in R/original_simulation.R
\name{simulate_seiqhrf}
\alias{simulate_seiqhrf}
\title{SEIQHRF Simulation Wrapper}
@@ -67,7 +67,8 @@ simulate_seiqhrf(
dq.rate = (7/365)/1000,
dh.rate = (20/365)/1000,
dr.rate = (7/365)/1000,
- out = "mean"
+ out = "mean",
+ seed = NULL
)
}
\arguments{
@@ -340,6 +341,8 @@ to model emigration as well as deaths.}
\item{out}{Summary function for the simulation runs. median is
also available, or percentiles, see the EpiModel documentation.}
+
+\item{seed}{Random see for checking consistency.}
}
\value{
list with simulation and simulation dataframe
diff --git a/man/summary_seiqhrf.Rd b/man/summary.seiqhrf.Rd
similarity index 75%
rename from man/summary_seiqhrf.Rd
rename to man/summary.seiqhrf.Rd
index 27b85cb..3e68728 100644
--- a/man/summary_seiqhrf.Rd
+++ b/man/summary.seiqhrf.Rd
@@ -1,13 +1,15 @@
% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/summary.R
-\name{summary_seiqhrf}
-\alias{summary_seiqhrf}
+% Please edit documentation in R/FIN_summary.R
+\name{summary.seiqhrf}
+\alias{summary.seiqhrf}
\title{summary function}
\usage{
-summary_seiqhrf(object)
+\method{summary}{seiqhrf}(object, ...)
}
\arguments{
-\item{object}{ICM parameters.}
+\item{object}{seiqhrf object, returned by \code{\link{seiqhrf}}.}
+
+\item{...}{Additional parameters.}
}
\value{
A list of compartments, each compartment is a list of three:
diff --git a/tests/testthat/.DS_Store b/tests/testthat/.DS_Store
index 41c1f5f62a267c92f511c48760e81d53b16b7820..9e1fc75f6bee687fada1a430d6e0d1018cdd04a1 100644
GIT binary patch
literal 6148
zcmeHKOHRWu5Pfcggg}DA29}XnuqqNah)P+og$vNHh}05DKN7IXJvatO;TRkO-q>Ea
zPS_#@Gu6!V?~R?;#;yUF-fDafbO3DA1(O|y29tjAnzek3bQ+tZivec*bTO}u)&^BT
z75Hrm$h&)pB_>!fhVO5IF;;w+rRNOSSF>U?oE1Zkm(AjYe)?&RC18IeM~M<S)^;~&
zVt+*Jy?{ODj8igY?qz*u^A>RrBJKmZCEHevjcaecL>W3kh8*{teS#6&aulr9p1pk$
z@xv7eKet@j33EI%S1m)b8)XR9lU0v6^O{)Vse%0-v5$iFFUZRAh&ytp{5}_~&8?OL
zJF2#pt=A|Ut7`v30iM}n+wyBvQw3B3RUlMA_J@=%n0c%m+O31ZUIB;+o2{`de-DI{
zddxgl4%tI9j!JY?<F6RT(b=B*ahb=;p`*k2%ZKqZ8-GJFc6R1Z?QodPp+;3e71&gu
z?W!F)|1ZD4|8J7CrV6M6e@X$9^l$n-K9WCMYY!)9ZA!nTi^;gkp`|eR<JcDDC_biJ
ZV_zx-V&<`O$QGLZ5wJ37Pz8QefiDKCkK+IU
delta 98
zcmZoMXfc=|#>CJ*u~2NHo}wrt0|NsP3otO~GGsF3F{Cl1G9*tdRA*$I%)@HL$KVH)
s*Z&U%n*}+<ST;7OGH+(*;O77u0u=bpJegm_kda}s4UhEZ7?Bms0K%^tfB*mh
diff --git a/tests/testthat/test-seiqhrf.R b/tests/testthat/test-all-modules.R
similarity index 91%
rename from tests/testthat/test-seiqhrf.R
rename to tests/testthat/test-all-modules.R
index 615b79d..acda5c0 100644
--- a/tests/testthat/test-seiqhrf.R
+++ b/tests/testthat/test-all-modules.R
@@ -1,4 +1,4 @@
-test_that("Identical output as Churches' original function: seiqhrf", {
+test_that("All modules(MOD) produce the same result as Churches' modules (.seiqhrf.icm)", {
full_params <- set_param(s.num = 1000, nsteps = 10, arec.rate = 0)
param <- full_params$param
diff --git a/tests/testthat/test-plot-model.R b/tests/testthat/test-plot-model.R
new file mode 100644
index 0000000..a42de96
--- /dev/null
+++ b/tests/testthat/test-plot-model.R
@@ -0,0 +1,49 @@
+test_that("plot_model is producing the same plot", {
+
+s.num <- 2000 # number susceptible
+i.num <- 15 # number infected
+q.num <- 5 # number in self-isolation
+h.num <- 1 #
+vals <- c(10, 7)
+timing <- c(7, 14)
+nsteps <- 90
+
+
+## baseline
+baseline_sim <- simulate_seiqhrf(nsteps = nsteps, s.num = s.num, i.num = i.num,
+ q.num = q.num, h.num = h.num)
+
+## Experiment 1
+act_rate <- vary_param(nstep = nsteps, vals = vals, timing = timing)
+
+sim_exp <- simulate_seiqhrf(nsteps = nsteps, s.num = s.num, i.num = i.num,
+ q.num = q.num, h.num = h.num, act.rate.e = act_rate,
+ act.rate.i = act_rate * 0.5)
+
+# Experiment #2
+vals <- c(10, 7, 7, 10)
+timing <- c(7, 14, 21, 28)
+act_rate_relax <- vary_param(nsteps, vals = vals, timing = timing)
+
+sim_exp_relax <- simulate_seiqhrf(nsteps = nsteps, s.num = s.num, i.num = i.num,
+ q.num = q.num, h.num = h.num,
+ act.rate.e = act_rate_relax,
+ act.rate.i = 0.5 * act_rate_relax)
+
+
+
+ori_plot <- plot_models(c(baseline_sim, sim_exp, sim_exp_relax),
+ sim_id = c('Baseline', 'Closures', 'Closures (2 mo)'),
+ start_date = lubridate::ymd("2020-01-01"),
+ comp_remove = c('s.num', 'r.num'),
+ plot_title = 'Closures Experiment')
+
+my_plot <- plot(list('Baseline' = baseline_sim$sim, 'Closures' = sim_exp$sim, 'Closures (2 mo)' = sim_exp_relax$sim),
+ start_date = lubridate::ymd("2020-01-01"),
+ comp_remove = c('s.num', 'r.num'),
+ plot_title = 'Closures Experiment')
+
+
+expect_equal(all.equal(ori_plot, my_plot), TRUE)
+
+})
diff --git a/tests/testthat/test-seiqhrf-simulate.R b/tests/testthat/test-seiqhrf-simulate.R
new file mode 100644
index 0000000..cb9158e
--- /dev/null
+++ b/tests/testthat/test-seiqhrf-simulate.R
@@ -0,0 +1,32 @@
+test_that("seiqhrf and simulate_seiqhrf produce identical output", {
+
+ s.num = 1000
+ q.num = 10
+ nsteps = 10
+ nsims = 3
+ arec.rate = 0 ### has to be fixed 0 for comparison
+
+ No_seeds <- 10
+ seed_list <- sample(1:1000, No_seeds)
+ comp <- rep(NA, No_seeds)
+ i <- 1
+ for(seed in seed_list){
+ Churhes_res <- simulate_seiqhrf(nsteps = nsteps, nsims = nsims,
+ arec.rate = arec.rate, s.num = s.num, q.num = q.num,
+ infection.FUN = infection.seiqhrf.icm,
+ recovery.FUN = progress.seiqhrf.icm,
+ departures.FUN = departures.seiqhrf.icm,
+ arrivals.FUN = arrivals.seiqhrf.icm, seed = seed)$sim
+ class(Churhes_res) <- "seiqhrf"
+
+ param <- param_seiqhrf(arec.rate = arec.rate)
+ init <- init_seiqhrf(s.num = s.num, q.num = q.num)
+ control <- control_seiqhrf(nsteps = nsteps, nsims = nsims)
+ sirplus_res <- seiqhrf(init, control, param, seed)
+
+ comp[i] <- identical(Churhes_res, sirplus_res)
+ i <- i + 1
+ }
+
+ expect_equal(sum(comp), No_seeds)
+})
diff --git a/vignettes/sirplus_intro.Rmd b/vignettes/sirplus_intro.Rmd
index 6dec74f..07d04fe 100644
--- a/vignettes/sirplus_intro.Rmd
+++ b/vignettes/sirplus_intro.Rmd
@@ -75,8 +75,7 @@ functions to examine the distributions of timings for various transitions based
on these parameters. In the case of a disease with observed data available, these plots can be used to sanity check parameter settings.
```{r baseline sims, fig.height=8}
-times <- get_times(baseline_sim)
-plot_times(times)
+plot_times(baseline_sim$sim)
```
diff --git a/vignettes/sirplus_intro_v2.Rmd b/vignettes/sirplus_intro_v2.Rmd
index 4a9bdd7..2594e00 100644
--- a/vignettes/sirplus_intro_v2.Rmd
+++ b/vignettes/sirplus_intro_v2.Rmd
@@ -16,15 +16,54 @@ knitr::opts_chunk$set(echo = TRUE, cache=FALSE, eval=TRUE, collapse = TRUE,
-```{r, set parameters}
+The sirplus package makes it easy to generate stochastic individual compartment
+models (ICMs) to simulate contagious disease spread using compartments not
+available in standard SIR packages. This method and most of the code was
+originally written by Tim Churches (see his [blog post](https://timchurches.github.io/blog/posts/2020-03-18-modelling-the-effects-of-public-health-interventions-on-covid-19-transmission-part-2/)). The sirplus package was developed by
+the Bioinformatics & Cellular Genomics team at St. Vincent's Institute of
+Medical Research in order to help St. Vincents' Hospital model the COVID-19
+pandemic.
+
+The compartments available in this package include:
+
+| Compartment | Functional definition |
+|-------------|-----------------------------------------------------------------------------------|
+| S | Susceptible individuals |
+| E | Exposed **and** infected, not yet symptomatic but potentially infectious |
+| I | Infected, symptomatic **and** infectious |
+| Q | Infectious, but (self-)isolated |
+| H | Requiring hospitalisation (would normally be hospitalised if capacity available) |
+| R | Recovered, immune from further infection |
+| F | Case fatality (death due to infection, not other causes) |
+
+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)
-s.num <- 1000 # number susceptible
+devtools::load_all()
+```
+
+
+## Simulate and inspect a baseline sirplus model
+
+### Set parameters
+
+Here we will simulate the epidemiological data for a made-up population with 1000
+susceptible individuals (S), 50 that are infected but not in the hospital or in
+self-quarantine (I; maybe people that are infected/symptomatic but not tested/
+aware), 10 confirmed cases that have self-isolated (Q), and 1 confirmed case that
+has been hospitalized (H). We call this the baseline model because it uses
+default parameters for disease spread (i.e. no additional interventions).
+
+```{r, set parameters}
+s.num <- 2000 # number susceptible
i.num <- 15 # number infected
q.num <- 5 # number in self-isolation
h.num <- 1 # number in the hospital
nsteps <- 90 # number of steps (e.g. days) to simulate
control <- control_seiqhrf(nsteps = nsteps)
-param <- param_seiqhrf(arec.rate = 0)
+param <- param_seiqhrf()
init <- init_seiqhrf(s.num = s.num, i.num = i.num, q.num = q.num, h.num = h.num)
print(init)
@@ -32,7 +71,42 @@ print(control)
print(param)
```
+### Simulate baseline
+
+This will produce an seiqhrf object.
+
```{r, run simulation}
sim <- seiqhrf(init, control, param)
sim
-```
\ No newline at end of file
+```
+
+### Extract summary of the simulation
+
+```{r, extract summary}
+res <- summary(sim)
+names(res)
+```
+
+### Inspect baseline transition distributions
+
+The sirplus model controls transitions between compartments, i.e. a change in state for an individual (e.g. going from self-isolation to hospital), using a variety of
+transition parameters. You can use the `plot()` functions, and set parameter `method` as `times`
+to examine the distributions of timings for various transitions based
+on these parameters. In the case of a disease with observed data available, these plots can be used to sanity check parameter settings.
+
+
+```{r baseline sims, fig.height=8}
+plot(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 by changing parameter `method` to `models`.
+
+```{r viz prevalence}
+plot(sim,
+ start_date = lubridate::ymd("2020-01-01"),
+ comp_remove = c('s.num', 'r.num'),
+ plot_title = 'Baseline Model')
+```
--
GitLab