diff --git a/.Rbuildignore b/.Rbuildignore
index 708610f47beec04f16f9dfe7f49791ded3b7bfb9..7a4f21b38a61a1c866165a1d84ff0d40d07f71a7 100644
--- a/.Rbuildignore
+++ b/.Rbuildignore
@@ -11,3 +11,4 @@
^pkgdown$
^doc$
^Meta$
+^\.vscode
diff --git a/DESCRIPTION b/DESCRIPTION
index 7de52d108555979a596b6c062a32d7b6760a0147..09f7f2e833156c5286238f1fce53cc35a8654339 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -18,27 +18,33 @@ Authors@R:
Description: sirplus is a package for the modeling of COVID-19 spread with extensions to the classical Susceptible-Infectious-Recovered (SIR) model, using
stochastic individual compartment models (ICMs). It provides a simple interface
for creating experiments to demonstrate how factors like social-distancing and
- epidemeological parameters will change the curve of an epidemic's trajectory.
+ epidemiological parameters will change the curve of an epidemic's trajectory.
License: GPL-3 + file LICENSE
LazyData: TRUE
Depends:
R (>= 3.6)
Imports:
- ggplot2,
- EpiModel,
- foreach,
doParallel,
dplyr,
tidyr,
stats,
future
+ EpiModel,
+ ggplot2,
+ foreach
Suggests:
+ earlyR,
+ knitr,
+ lubridate,
magick,
- testthat,
rmarkdown,
spelling,
+<<<<<<< HEAD
earlyR,
knitr
+=======
+ testthat
+>>>>>>> 382b87aca1de22f5aeea1065b2d6138cc69bc13a
biocViews: Epidemiology, Software
RoxygenNote: 7.1.0
VignetteBuilder: knitr
diff --git a/NAMESPACE b/NAMESPACE
index 9b3b2eda5f0ba94bb00d0352e8ffa2116b03928d..17281a0820400b94d0be14ab47492b48d1aa3673 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -37,3 +37,4 @@ importFrom(stats,rbinom)
importFrom(stats,rgeom)
importFrom(stats,sd)
importFrom(tidyr,pivot_longer)
+importFrom(tidyr,pivot_wider)
diff --git a/R/get_weekly_local.R b/R/get_weekly_local.R
index b98a9425cfe929e95e992f181eabd339cdbcfe38..77ab396abc54e086c9ca083cef5cd76c9499e7bc 100644
--- a/R/get_weekly_local.R
+++ b/R/get_weekly_local.R
@@ -2,58 +2,89 @@
#'
#'
#' @param sim An \code{icm} 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 hospitality
-#' @param total_population True population size, needed only if simulation size is smaller than the true population size due to computational cost etc.
+#' @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 hospitality and ICU respectively
+#' \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 hospitality locally,
-#' \item \code{icu:} the number of patients that require ICU locally. }
+#' \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. }
#
#' }
#'
-
-
-
-get_weekly_local <- function(sim, market.share = .4, icu_percent = .1, total_population = NULL){
+#' @importFrom tidyr pivot_wider
+#'
+get_weekly_local <- function(sim,
+ market.share = .04,
+ icu_percent = .1,
+ start_date = 'na',
+ time_lim = 90,
+ total_population = NULL){
hosp <- sim$df$h.num
-
if(!is.null(total_population)){
- if(total_population < max(sim$df$s.num)) stop("total Population should be larger than simulated size")
+ if(total_population < max(sim$df$s.num))
+ stop("total Population should be larger than simulated size")
cat("Scalling w.r.t total population")
hosp <- hosp*total_population/max(sim$df$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")
+ 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")
hosp[is.na(hosp)] <- 0
+ hosp <- hosp[1: time_lim]
hosp_week <- split(hosp, ceiling(seq_along(hosp)/7))
hosp_sum_week <- unlist(lapply(hosp_week, sum))
t_sz <- length(hosp_sum_week)
- plot_hosp_icu_week <- data.frame(wk = rep(seq_along(hosp_sum_week), 2),
- hosp_icu = c(hosp_sum_week, hosp_sum_week*icu_percent),
- group = rep(c("general hopitality", "icu"), each = t_sz))
-
+ hosp_wk_df <- data.frame(wk = rep(seq_along(hosp_sum_week), 2),
+ group = rep(c("general", "icu"),
+ each = t_sz),
+ hosp_icu = c(hosp_sum_week -
+ (hosp_sum_week*icu_percent),
+ hosp_sum_week*icu_percent))
- gg <- ggplot(data=plot_hosp_icu_week, aes(x = wk, y = hosp_icu, fill = group)) +
- geom_bar(stat="identity") +
- labs(y="Counts", x = "Week") +
- scale_x_continuous(breaks = seq(0,t_sz,5), labels= seq(0,t_sz,5))
+ if(class(start_date) == 'Date'){
+
+ hosp_wk_df <- data.frame(append(hosp_wk_df,
+ list(Date=start_date +
+ (7 * (hosp_wk_df$wk - 1))),
+ after=match("wk", names(hosp_wk_df))))
+
+ gg <- ggplot(data=hosp_wk_df, aes(x = Date, y = hosp_icu, fill = group)) +
+ geom_bar(stat="identity") + theme_bw() +
+ scale_x_date(date_breaks = "1 week", date_labels = "%m-%d") +
+ labs(y="Weekly Hospital Load (sum over week)", x = "Week")
+
+ }else{
+
+ gg <- ggplot(data=hosp_wk_df, aes(x = wk, y = hosp_icu, fill = group)) +
+ geom_bar(stat="identity") + theme_bw() +
+ labs(y="Weekly Hospital Load (sum over week)", x = "Week") +
+ scale_x_continuous(breaks = seq(0,t_sz,5), labels= seq(0,t_sz,5))
+ }
- res <- data.frame(wk = seq_along(hosp_sum_week), hosp = hosp_sum_week, icu = hosp_sum_week*icu_percent)
- rownames(res) <- c()
+ res <- hosp_wk_df %>% tidyr::pivot_wider(names_from = group, values_from = hosp_icu)
return(list("plot" = gg, "result" = res))
-
+
}
diff --git a/R/results-parse.R b/R/results-parse.R
index 1ca51148c94db71084db607fcc9e6a39169c6e64..843e6c9a38c7535c84a6ead62c509398a3a87eb4 100644
--- a/R/results-parse.R
+++ b/R/results-parse.R
@@ -150,6 +150,7 @@ plot_models <- function(sims = baseline_sim,
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()
@@ -157,6 +158,7 @@ plot_models <- function(sims = baseline_sim,
plot_df %>% ggplot(aes(x = Date, y = count, colour = compartment,
linetype = sim)) +
facet_grid(reo_exp(experiment) ~ .) +
+ 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") +
@@ -167,6 +169,7 @@ plot_models <- function(sims = baseline_sim,
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") +
@@ -176,6 +179,7 @@ plot_models <- function(sims = baseline_sim,
linetype = sim)) +
facet_grid(reo_exp(experiment) ~ .) +
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") +
diff --git a/README.md b/README.md
index 880f2ce3d5abdef1821b8b62d390c60b9379ecdc..8cf2de04a752f71fe612e83df0bbf76be72d30e4 100644
--- a/README.md
+++ b/README.md
@@ -1,7 +1,6 @@
# sirplus
-sirplus is a package for the modeling of COVID-19 spread using stochastic individual compartment models (ICMs). The model implemented extends the classical Susceptible-Infectious-Recovered (SIR) model by adding compartments for Exposed, Quarantined, Hospitalised and Case Fatality individuals. The package provides a simple interface for creating experiments to demonstrate how factors like social-distancing and epidemeological parameters will change the curve.
-
+sirplus is a package for the modeling of COVID-19 spread using stochastic individual compartment models (ICMs). The model implemented extends the classical Susceptible-Infectious-Recovered (SIR) model by adding compartments for Exposed, Quarantined, Hospitalised and Case Fatality individuals. The package provides a simple interface for creating experiments to demonstrate how factors like social-distancing and epidemiological parameters will change the curve.
To view the vignette install the package then:
diff --git a/man/get_weekly_local.Rd b/man/get_weekly_local.Rd
index a906e9e3949a249aff826091a1d8da811aea09e2..abfb0880a6d2bd5501f032d5ed0495b7e166d388 100644
--- a/man/get_weekly_local.Rd
+++ b/man/get_weekly_local.Rd
@@ -6,28 +6,40 @@
\usage{
get_weekly_local(
sim,
- market.share = 0.4,
+ market.share = 0.04,
icu_percent = 0.1,
+ start_date = "na",
+ time_lim = 90,
total_population = NULL
)
}
\arguments{
\item{sim}{An \code{icm} object returned by \link{simulate_seiqhrf}.}
-\item{market.share}{between 0 and 1, percentage of local hospital beds in the simulated unit (e.g. state)}
+\item{market.share}{between 0 and 1, percentage of local hospital beds in
+the simulated unit (e.g. state)}
-\item{icu_percent}{between 0 and 1, percentage of patients that should go to ICU among the ones that need hospitality}
+\item{icu_percent}{between 0 and 1, percentage of patients that should go to
+ICU among the ones that need hospitalization}
-\item{total_population}{True population size, needed only if simulation size is smaller than the true population size due to computational cost etc.}
+\item{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.}
+
+\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{
-\item \code{plot:} A \code{ggplot} object, bar charts of count of patients requiring hospitality and ICU respectively
+\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 hospitality locally,
- \item \code{icu:} the number of patients that require ICU locally. }
+ \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. }
}
}
\description{
diff --git a/man/simulate.seiqhrf.Rd b/man/simulate.seiqhrf.Rd
deleted file mode 100644
index 824e0c9b0317c857e79a3b64e1a03baa475e2310..0000000000000000000000000000000000000000
--- a/man/simulate.seiqhrf.Rd
+++ /dev/null
@@ -1,299 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/simulation.R
-\name{simulate.seiqhrf}
-\alias{simulate.seiqhrf}
-\title{SEIQHRF Simulation Wrapper}
-\usage{
-\method{simulate}{seiqhrf}(
- type = "SEIQHRF",
- nsteps = 366,
- nsims = 8,
- ncores = 4,
- prog.rand = FALSE,
- rec.rand = FALSE,
- fat.rand = TRUE,
- quar.rand = FALSE,
- hosp.rand = FALSE,
- disch.rand = TRUE,
- infection.FUN = "infection.FUN",
- recovery.FUN = "recovery.FUN",
- departures.FUN = "departures.FUN",
- arrivals.FUN = "arrivals.FUN",
- get_prev.FUN = "get_prev.FUN",
- s.num = 9997,
- e.num = 0,
- i.num = 3,
- q.num = 0,
- h.num = 0,
- r.num = 0,
- f.num = 0,
- inf.prob.e = 0.02,
- act.rate.e = 10,
- inf.prob.i = 0.05,
- act.rate.i = 10,
- inf.prob.q = 0.02,
- act.rate.q = 2.5,
- quar.rate = 1/30,
- hosp.rate = 1/100,
- disch.rate = 1/15,
- prog.rate = 1/10,
- prog.dist.scale = 5,
- prog.dist.shape = 1.5,
- rec.rate = 0.071,
- rec.dist.scale = 35,
- rec.dist.shape = 1.5,
- fat.rate.base = 1/50,
- hosp.cap = 40,
- fat.rate.overcap = 1/25,
- fat.tcoeff = 0.5,
- vital = TRUE,
- a.rate = (10.5/365)/1000,
- a.prop.e = 0.01,
- a.prop.i = 0.001,
- a.prop.q = 0.01,
- ds.rate = (7/365)/1000,
- de.rate = (7/365)/1000,
- di.rate = (7/365)/1000,
- dq.rate = (7/365)/1000,
- dh.rate = (20/365)/1000,
- dr.rate = (7/365)/1000,
- out = "mean"
-)
-}
-\arguments{
-\item{type}{Type of model: SI, SIR, SIS, SEIR, SEIQHR and
-SEIQHRF available, but only SEIQHRF is likely to work in the
-current version of the code.}
-
-\item{nsteps}{Number of time steps to solve the model over. This must be a
-positive integer.}
-
-\item{nsims}{Number of simulations to run.}
-
-\item{ncores}{Number of physical CPU cores used for parallel computation.}
-
-\item{prog.rand}{Method for progression from E compartment to I. If TRUE,
-random binomial draws at prog.rate, if FALSE, random draws from a
-Weibull distribution (yes, I know it should be a discrete Weibull
-distribution but it makes little difference and speed of computation
-matters), with parameters prog.dist.scale and prog.dist.shape}
-
-\item{rec.rand}{Method for recovery transition from I, Q or H to R. If TRUE,
-random binomial draws at prog.rate, if FALSE, random draws from a
-random draws from a Weibull distribution, with parameters
-rec.dist.scale and rec.dist.shape.}
-
-\item{fat.rand}{Method for case fatality transition from H to F. If TRUE,
-random binomial draws at fat.rate.base, if FALSE, random sample with
-a sample fraction also given by fat.rate.base. However, if the
-current number of patients in the H (needs hospitalisation)
-compartment is above a hospital capacity level specified by
-hosp.cap, then the fatality rate is the mean of the base fatality
-rate weighted by the hospital capacity, plus a higher rate,
-specified by fat.rate.overcap, weighted by the balance of those
-requiring hospitalisation (but who can't be accommodated). By
-setting fat.rate.overcap higher, the effect of exceeding the
-capacity of the health care system can be simulated. There is also
-a coefficient fat.tcoeff for the fatality rates that increases them
-as a linear function of the number of days the individual has been
-in the H compartment. Use of the co-efficient better approximates
-the trapezoid survival time distribution typical of ICU patients.}
-
-\item{quar.rand}{Method for self-isolation transition from I to Q. If TRUE,
-random binomial draws at quar.rate, if FALSE, random sample with a
-sample fraction also given by `quar.rate.}
-
-\item{hosp.rand}{Method for transition from I or Q to H -- that is, from
-infectious or from self-isolated to requiring hospitalisation. If
-TRUE, random binomial draws at hosp.rate, if FALSE, random sample
-with a sample fraction also given by `hosp.rate.}
-
-\item{disch.rand}{Method for transition from H to R -- that is, from
-requiring hospitalisation to recovered. If TRUE, random binomial
-draws at disch.rate, if FALSE, random sample with a sample fraction
-also given by disch.rate. Note that the only way out of the H
-compartment is recovery or death.}
-
-\item{infection.FUN}{No, being infected with SARS-CoV2 is not fun. Rather
-this is the the name of the function to implement infection
-processes. Use the default.}
-
-\item{recovery.FUN}{Function to implement recovery processes. Use the
-default.}
-
-\item{departures.FUN}{Handles background demographics, specifically
-departures (deaths not due to the virus, and emigration). Use the
-default.}
-
-\item{arrivals.FUN}{Handles background demographics, specifically arrivals
-(births and immigration). Uses the original EpiModel code
-currently. A replacement that implements modelling the arrival of
-infected individuals is under development -- but for now, all
-arrivals go into the S compartment.}
-
-\item{get_prev.FUN}{Utility function that collects prevalence and transition
-time data from each run and stores it away in the simulation result
-object. Use the default.}
-
-\item{s.num}{Initial number of *S compartment individuals in
-the simulated population. An overall population of 10,000 is a good
-compromise. A set of models will still take several minutes or more
-to run, in parallel.}
-
-\item{e.num}{Initial number of E compartment individuals in
-the simulated population.}
-
-\item{i.num}{Initial number of I compartment individuals in
-the simulated population.}
-
-\item{q.num}{Initial number of Q compartment individuals in
-the simulated population.}
-
-\item{h.num}{Initial number of H compartment individuals in
-the simulated population.}
-
-\item{r.num}{Initial number of R compartment individuals in
-the simulated population.}
-
-\item{f.num}{Initial number of F compartment individuals in
-the simulated population.}
-
-\item{inf.prob.e}{Probability of passing on infection at each
-exposure event for interactions between infectious people in the E
-compartment and susceptibles in S. Note the default is lower than for
-inf.prob.i reflecting the reduced infectivity of infected but
-asymptomatic people (the E compartment). Otherwise as for inf.exp.i.}
-
-\item{act.rate.e}{The number of exposure events (acts) between
-infectious individuals in the E compartment and susceptible individuals
-in the S compartment, per day. Otherwise as for act.rate.i.}
-
-\item{inf.prob.i}{Probability of passing on infection at each
-exposure event for interactions between infectious people in the I
-compartment and susceptibles in S. Reducing inf.prob.i is equivalent to
-increasing hygiene measures, such as not putting hands in eyes, nose or
-moth, use of hand sanitisers, wearing masks by the infected, and so on.}
-
-\item{act.rate.i}{The number of exposure events (acts) between
-infectious individuals in the I compartment and susceptible individuals
-in the S compartment, per day. It's stochastic, so the rate is an
-average, some individuals may have more or less. Note that not every
-exposure event results in infection - that is governed by the inf.prob.i
-parameters (see below). Reducing act.rate.i is equivalent to increasing
-social distancing by people in the I compartment.}
-
-\item{inf.prob.q}{Probability of passing on infection at each
-exposure event for interactions between infectious people in the Q
-compartment and susceptibles in S. Note the default is lower than for
-inf.prob.i reflecting the greater care that self-isolated individuals
-will, on average, take regarding hygiene measures, such as wearing masks,
-to limit spread to others. Otherwise as for inf.exp.i.}
-
-\item{act.rate.q}{The number of exposure events (acts) between
-infectious individuals in the Q compartment (isolated, self or otherwise)
-and susceptible individuals in the S compartment, per day. Note the much
-lower rate than for the I and E compartments, reflecting the much
-greater degree of social isolation for someone in (self-)isolation. The
-exposure event rate is not zero for this group, just much less.
-Otherwise as for act.rate.i.}
-
-\item{quar.rate}{Rate per day at which symptomatic (or tested
-positive), infected I compartment people enter self-isolation (Q
-compartment). Asymptomatic E compartment people can't enter
-self-isolation because they don't yet know they are infected. Default is
-a low rate reflecting low community awareness or compliance with
-self-isolation requirements or practices, but this can be tweaked when
-exploring scenarios.}
-
-\item{hosp.rate}{Rate per day at which symptomatic (or tested
-positive), infected I compartment people or self-isolated Q compartment
-people enter the state of requiring hospital care -- that is, become
-serious cases. A default rate of 1% per day with an average illness
-duration of about 10 days means a bit less than 10% of cases will
-require hospitalisation, which seems about right (but can be tweaked,
-of course).}
-
-\item{disch.rate}{Rate per day at which people needing
-hospitalisation recover.}
-
-\item{prog.rate}{Rate per day at which people who are infected
-but asymptomatic (E compartment) progress to becoming symptomatic (or
-test-positive), the I compartment. See prog.rand above for more details.}
-
-\item{prog.dist.scale}{Scale parameter for Weibull distribution
-for progression, see prog.rand for details.}
-
-\item{prog.dist.shape}{Shape parameter for Weibull distribution
-for progression, see prog.rand for details. Read up on the Weibull
-distribution before changing the default.}
-
-\item{rec.rate}{Rate per day at which people who are infected and
-symptomatic (I compartment) recover, thus entering the R compartment.
-See rec.rand above for more details.}
-
-\item{rec.dist.scale}{Scale parameter for Weibull distribution for
-recovery, see rec.rand for details.}
-
-\item{rec.dist.shape}{Shape parameter for Weibull distribution for
-recovery, see rec.rand for details. Read up on the Weibull distribution
-before changing the default.}
-
-\item{fat.rate.base}{Baseline mortality rate per day for people
-needing hospitalisation (deaths due to the virus). See fat.rand for more
-details.}
-
-\item{hosp.cap}{Number of available hospital beds for the modelled
-population. See fat.rand for more details.}
-
-\item{fat.rate.overcap}{Mortality rate per day for people needing
-hospitalisation but who can't get into hospital due to the hospitals
-being full (see hosp.cap and fat.rand). The default rate is twice that
-for those who do get into hospital.}
-
-\item{fat.tcoeff}{Time co-efficient for increasing mortality rate
-as time in the H compartment increases for each individual in it. See
-fat.rand for details.}
-
-\item{vital}{Enables demographics, that is, arrivals and
-departures, to and from the simulated population.}
-
-\item{a.rate}{Background demographic arrival rate. Currently all
-arrivals go into the S compartment, the default is approximately the
-daily birth rate for Australia. Will be extended to cover immigration in
-future versions.}
-
-\item{a.prop.e}{Arrivals proportion that goes to E (immigration).}
-
-\item{a.prop.i}{Arrivals proportion that goes to I (immigration).}
-
-\item{a.prop.q}{Arrivals proportion that goes to Q (immigration).}
-
-\item{ds.rate}{Background demographic departure (death not due to
-virus) rates. Defaults based on Australian crude death rates. Can be
-used to model emigration as well as deaths.}
-
-\item{de.rate}{Background demographic departure (death not due to
-virus) rates. Defaults based on Australian crude death rates. Can be used
-to model emigration as well as deaths.}
-
-\item{dq.rate}{Background demographic departure (death not due to
-virus) rates. Defaults based on Australian crude death rates. Can be used
-to model emigration as well as deaths.}
-
-\item{dh.rate}{Background demographic departure (death not due to
-virus) rates. Defaults based on Australian crude death rates. Can be used
-to model emigration as well as deaths.}
-
-\item{dr.rate}{Background demographic departure (death not due to
-virus) rates. Defaults based on Australian crude death rates. Can be used
-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.}
-}
-\value{
-list with simulation and simulation dataframe
-}
-\description{
-Wrapper to implement SEIQHRF model
-}
diff --git a/vignettes/sirplus_intro.Rmd b/vignettes/sirplus_intro.Rmd
index 6592f47e2b1473745923ec4e69c0d47c517da223..150233f3288cb832d2cdfada79f536d4e958cdce 100644
--- a/vignettes/sirplus_intro.Rmd
+++ b/vignettes/sirplus_intro.Rmd
@@ -1,8 +1,8 @@
---
title: "Introduction to sirplus"
-date: "Last updated: 23 March 2020"
+date: "Last updated: 26 March 2020"
vignette: >
- %\VignetteIndexEntry{introduction to sirplus}
+ %\VignetteIndexEntry{Introduction to sirplus}
%\VignetteEngine{knitr::rmarkdown}
\usepackage[utf8]{inputenc}
---
@@ -13,15 +13,15 @@ knitr::opts_chunk$set(echo = TRUE, cache=FALSE, eval=TRUE, collapse = TRUE,
tidy.opts=list(width.cutoff=60), tidy=FALSE)
```
-
+
The sirplus package makes it easy to generate stochastic individual compartment
-models (ICMs) to simulate contageous disease spread using compartments not
+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. Vincents' Institute of
+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.
+pandemic.
The compartments available in this package include:
@@ -30,29 +30,28 @@ The compartments available in this package include:
| 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) |
+| 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 COVID-19, not other causes) |
+| 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 [blog post](https://timchurches.github.io/blog/posts/2020-03-18-modelling-the-effects-of-public-health-interventions-on-covid-19-transmission-part-2/).
-
+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)
devtools::load_all(".")
```
-## Simulate and Inspect a Baseline sirplus Model
+## Simulate and inspect a baseline sirplus model
-### Simulate
+### Simulate
-Here we will simulate the epi data for a made up population with 1000
-susceptible individuals, 50 that are infected but not in the hospital or in
-self-quarentine (maybe people that are infected/symptomatic but not tested/
-aware), 10 confirmed cases that have self-isolated, and 1 confirmed case that
-has been hospitalized. We call this the baseline model because it uses
+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 simulate baselines}
@@ -61,30 +60,29 @@ i.num <- 50 # number infected
q.num <- 10 # number in self-isolation
h.num <- 1 # number in the hospital
-baseline_sim <- simulate_seiqhrf(s.num = s.num, i.num = i.num,
+baseline_sim <- simulate_seiqhrf(s.num = s.num, i.num = i.num,
q.num = q.num, h.num = h.num)
head(baseline_sim$df)
```
+### Inspect baseline transition distributions
-### Inspect Baseline Transition Distributions
-
-The sirplus model controls transitions between compartments using a variety of
+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 `get_times()` and `plot_times()`
functions to examine the distributions of timings for various transitions based
-on these parameters.
+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}
+```{r baseline sims, fig.height=8}
times <- get_times(baseline_sim)
plot_times(times)
```
-### Plot Baseline sirplus Results
+### Plot baseline sirplus results
-To visualise your sirplus model, 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.
```{r viz prevalence}
plot_models(sims = baseline_sim,
@@ -94,23 +92,17 @@ plot_models(sims = baseline_sim,
plot_title = 'Baseline Model')
```
-
## Run an experiment
-With the sirplus package you can also set up experiments. For example, lets say
-that one week after the beginning of the epidemic, schools and non-essential
-buisnesses are closed to encourage social distancing. We can model the impact
-that policy will have and compare it to our baseline by ramping down the
-act.rate.e (# of exposure events, or acts, between individuals in the E and S
-compartment, per day) after 7 days and plot this model next to our baseline.
-
+With the sirplus package you can also set up experiments. For example, let's say that one week after the beginning of the epidemic, schools and non-essential businesses are closed to encourage social distancing. We can model the impacts that policies may have and compare to our baseline model. In this example, we model a policy-based increase in social distancing by ramping down the act.rate.e (# of exposure events, or acts, between individuals in the E and S compartment, per day) after 7 days and plot this model next to our baseline.
```{r experiment 1}
closures_RampOnday7 <- function(t) {
ifelse(t < 7, 10, ifelse(t <= 14, 10 - (t-7)*(10 - 5)/7, 5))
}
-closures_sim <- simulate_seiqhrf(s.num = s.num, i.num = i.num,
+
+closures_sim <- simulate_seiqhrf(s.num = s.num, i.num = i.num,
q.num = q.num, h.num = h.num,
act.rate.e = closures_RampOnday7(1:366),
act.rate.i = closures_RampOnday7(1:366))
@@ -125,6 +117,28 @@ plot_models(sims = c(baseline_sim, closures_sim),
```
From these results we see that this policy would likey reduce the peak number
-of infections from 700 to 500 and would flatten the curve for infections and
-thus hospitalizations.
-
+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 models to aid hospitals
+
+Use the `get_weekly_local()` function to extract and plot the number of weekly expected
+patients in a hospital. This function takes the following input:
+
+ - `sims`: Single or list of outputs from simulate.seiqhrf().
+ - `market.share`: Percent of cases in your model population (s.num) that are
+ anticipated at the hospital of interest. Default = 4% (0.04).
+ - `icu_percent`: Percent of hospitalised cases are are likely to need treatment in an intensive care unit (ICU). Default = 10% (0.1).
+ - `start_date`: Epidemic start date. Default is `NA`.
+ - `time_limit`: Number of days (nsteps) to include. Default = 90.
+ - `total_population`: True population size. This parameter is only needed if
+ the simulation size (s.num) was smaller than the true population size (i.e.
+ scaled down) to reduce computational cost.
+
+```{r hospital visualization}
+weekly_local <- get_weekly_local(baseline_sim, time_lim = 50, start_date = lubridate::ymd("2020-01-01"))
+
+head(weekly_local$result)
+weekly_local$plot
+```