From 65058422de43212e493288a4775bcac831cd5e8e Mon Sep 17 00:00:00 2001
From: cazodi <cazodi@svi.edu.au>
Date: Tue, 24 Mar 2020 12:29:23 +1100
Subject: [PATCH] update description and namespace
---
DESCRIPTION | 98 +-
NAMESPACE | 175 ----
R/data.R | 76 --
R/epi-functions.R | 1524 +++++++++++++++++++++++++++++++
R/results-parse.R | 40 +
R/sirPlus-functions.R | 956 +++++++++++++++++++
R/sirPlus-package.R | 9 +
R/utils.R | 84 +-
man/arrivals.FUN.Rd | 23 +
man/control.icm.Rd | 39 +
man/cum_discr_si.Rd | 21 +
man/departures.FUN.Rd | 20 +
man/get_prev.FUN.Rd | 20 +
man/get_times.Rd | 17 +
man/icm.seiqhrf.Rd | 21 +
man/infection.FUN.Rd | 19 +
man/infection.seiqhrf.icm.Rd | 19 +
man/init_status.icm.Rd | 17 +
man/initialize.FUN.Rd | 21 +
man/merge.seiqhrf.icm.Rd | 19 +
man/progress.seiqhrf.icm.Rd | 19 +
man/saveout.seiqhrf.icm.Rd | 21 +
man/set.control.Rd | 54 ++
man/set.init.Rd | 47 +
man/set.param.Rd | 175 ++++
man/sirPlus.Rd | 10 +
vignettes/epimodels-SEIQHRF.Rmd | 133 +++
27 files changed, 3300 insertions(+), 377 deletions(-)
create mode 100644 R/epi-functions.R
create mode 100644 R/results-parse.R
create mode 100644 R/sirPlus-functions.R
create mode 100644 R/sirPlus-package.R
create mode 100644 man/arrivals.FUN.Rd
create mode 100644 man/control.icm.Rd
create mode 100644 man/cum_discr_si.Rd
create mode 100644 man/departures.FUN.Rd
create mode 100644 man/get_prev.FUN.Rd
create mode 100644 man/get_times.Rd
create mode 100644 man/icm.seiqhrf.Rd
create mode 100644 man/infection.FUN.Rd
create mode 100644 man/infection.seiqhrf.icm.Rd
create mode 100644 man/init_status.icm.Rd
create mode 100644 man/initialize.FUN.Rd
create mode 100644 man/merge.seiqhrf.icm.Rd
create mode 100644 man/progress.seiqhrf.icm.Rd
create mode 100644 man/saveout.seiqhrf.icm.Rd
create mode 100644 man/set.control.Rd
create mode 100644 man/set.init.Rd
create mode 100644 man/set.param.Rd
create mode 100644 man/sirPlus.Rd
create mode 100644 vignettes/epimodels-SEIQHRF.Rmd
diff --git a/DESCRIPTION b/DESCRIPTION
index d252d16..f8c1412 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,76 +1,52 @@
-Package: splatter
+Package: sirPLUS
Type: Package
-Title: Simple Simulation of Single-cell RNA Sequencing Data
-Version: 1.11.1
-Date: 2020-01-30
+Title: Using stochastic individual compartment models (ICMs) to simulate COVID-19 spread
+Version: 1.0
+Date: 2020-03-23
Authors@R:
- c(person("Luke", "Zappia", role = c("aut", "cre"),
- email = "luke@lazappi.id.au",
- comment = c(ORCID = "0000-0001-7744-8565")),
- person("Belinda", "Phipson", role = c("aut"),
- email = "belinda.phipson@mcri.edu.au",
- comment = c(ORCID = "0000-0002-1711-7454")),
- person("Alicia", "Oshlack", role = c("aut"),
- email = "alicia.oshlack@mcri.edu.au",
- comment = c(ORCID = "0000-0001-9788-5690")))
-Description: Splatter is a package for the simulation of single-cell RNA
- sequencing count data. It provides a simple interface for creating complex
- simulations that are reproducible and well-documented. Parameters can be
- estimated from real data and functions are provided for comparing real and
- simulated datasets.
+ c(person("Christina", "Azodi", role = c("aut"),
+ email = "cazodi@svi.edu.au",
+ comment = c(ORCID = "0000-0002-6097-606X")),
+ person("Puxue", "Qiao", role = c("aut"),
+ email = "pqiao@svi.edu.au"),
+ person("Ruqian", "Lyu", role = c("aut"),
+ email = "rlyu@svi.edu.au"),
+ person("Davis", "McCarthy", role = c("aut", "cre"),
+ email = "dmccarthy@svi.edu.au",
+ comment = c(ORCID = "0000-0002-2218-6833")))
+Description: sirPlus is a package for the modeling of COVID-19 spread 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.
License: GPL-3 + file LICENSE
LazyData: TRUE
Depends:
- R (>= 3.6),
- SingleCellExperiment
+ R (>= 3.6)
Imports:
- akima,
- BiocGenerics,
- BiocParallel,
- checkmate,
- edgeR,
- dplyr,
- data.table,
- fitdistrplus,
+ tidyverse,
+ magrittr,
+ lubridate,
+ stringr,
+ tibble,
+ broom,
ggplot2,
- locfit,
- matrixStats,
- methods,
- scales,
- scater,
- stats,
- SummarizedExperiment,
- utils,
- crayon,
- S4Vectors,
- rlang
-Suggests:
- BiocStyle,
- covr,
- cowplot,
+ gt,
knitr,
+ devtools,
+ DiagrammeR,
+ parallel,
+ foreach,
+ tictoc,
+ EpiModel,
+ incidence,
+ earlyR
+Suggests:
magick,
- limSolve,
- lme4,
- progress,
- pscl,
testthat,
rmarkdown,
- scDD,
- scran,
- mfa,
- phenopath,
- BASiCS (>= 1.7.10),
- zinbwave,
- SparseDC,
BiocManager,
- spelling,
- igraph,
- DropletUtils
-biocViews: SingleCell, RNASeq, Transcriptomics, GeneExpression, Sequencing,
- Software, ImmunoOncology
-URL: https://github.com/Oshlack/splatter
-BugReports: https://github.com/Oshlack/splatter/issues
+ spelling
+biocViews: Epidemiology, Software
RoxygenNote: 7.0.2
VignetteBuilder: knitr
Encoding: UTF-8
diff --git a/NAMESPACE b/NAMESPACE
index 758ebc8..1feab89 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -1,178 +1,3 @@
# Generated by roxygen2: do not edit by hand
-S3method(BASiCSEstimate,SingleCellExperiment)
-S3method(BASiCSEstimate,matrix)
-S3method(kersplatEstimate,SingleCellExperiment)
-S3method(kersplatEstimate,matrix)
-S3method(lun2Estimate,SingleCellExperiment)
-S3method(lun2Estimate,matrix)
-S3method(lunEstimate,SingleCellExperiment)
-S3method(lunEstimate,matrix)
-S3method(mfaEstimate,SingleCellExperiment)
-S3method(mfaEstimate,matrix)
-S3method(phenoEstimate,SingleCellExperiment)
-S3method(phenoEstimate,matrix)
-S3method(scDDEstimate,SingleCellExperiment)
-S3method(scDDEstimate,default)
-S3method(scDDEstimate,matrix)
-S3method(simpleEstimate,SingleCellExperiment)
-S3method(simpleEstimate,matrix)
-S3method(sparseDCEstimate,SingleCellExperiment)
-S3method(sparseDCEstimate,matrix)
-S3method(splatEstimate,SingleCellExperiment)
-S3method(splatEstimate,matrix)
-S3method(zinbEstimate,SingleCellExperiment)
-S3method(zinbEstimate,matrix)
-export(BASiCSEstimate)
-export(BASiCSSimulate)
-export(addGeneLengths)
-export(compareSCEs)
-export(diffSCEs)
-export(eQTLEstimate)
-export(eQTLSimulate)
-export(getParam)
-export(getParams)
-export(kersplatEstimate)
-export(kersplatSample)
-export(kersplatSetup)
-export(kersplatSimulate)
-export(listSims)
-export(lun2Estimate)
-export(lun2Simulate)
-export(lunEstimate)
-export(lunSimulate)
-export(makeCompPanel)
-export(makeDiffPanel)
-export(makeOverallPanel)
-export(mfaEstimate)
-export(mfaSimulate)
-export(newBASiCSParams)
-export(newKersplatParams)
-export(newLun2Params)
-export(newLunParams)
-export(newMFAParams)
-export(newPhenoParams)
-export(newSCDDParams)
-export(newSimpleParams)
-export(newSparseDCParams)
-export(newSplatParams)
-export(newZINBParams)
-export(neweQTLParams)
-export(phenoEstimate)
-export(phenoSimulate)
-export(scDDEstimate)
-export(scDDSimulate)
-export(setParam)
-export(setParams)
-export(simpleEstimate)
-export(simpleSimulate)
-export(sparseDCEstimate)
-export(sparseDCSimulate)
-export(splatEstimate)
-export(splatSimulate)
-export(splatSimulateGroups)
-export(splatSimulatePaths)
-export(splatSimulateSingle)
-export(splatSimulateeQTL)
-export(summariseDiff)
-export(zinbEstimate)
-export(zinbSimulate)
-exportClasses(BASiCSParams)
-exportClasses(KersplatParams)
-exportClasses(Lun2Params)
-exportClasses(LunParams)
-exportClasses(MFAParams)
-exportClasses(PhenoParams)
-exportClasses(SCDDParams)
-exportClasses(SimpleParams)
-exportClasses(SparseDCParams)
-exportClasses(SplatParams)
-exportClasses(ZINBParams)
-exportClasses(eQTLParams)
-importFrom(BiocParallel,SerialParam)
-importFrom(BiocParallel,bplapply)
-importFrom(S4Vectors,"metadata<-")
-importFrom(S4Vectors,metadata)
-importFrom(SingleCellExperiment,"cpm<-")
-importFrom(SingleCellExperiment,SingleCellExperiment)
-importFrom(SingleCellExperiment,cbind)
-importFrom(SingleCellExperiment,cpm)
-importFrom(SummarizedExperiment,"assays<-")
-importFrom(SummarizedExperiment,"colData<-")
-importFrom(SummarizedExperiment,"rowData<-")
-importFrom(SummarizedExperiment,assays)
-importFrom(SummarizedExperiment,colData)
-importFrom(SummarizedExperiment,rowData)
-importFrom(checkmate,checkCharacter)
-importFrom(checkmate,checkClass)
-importFrom(checkmate,checkDataFrame)
-importFrom(checkmate,checkFactor)
-importFrom(checkmate,checkFlag)
-importFrom(checkmate,checkInt)
-importFrom(checkmate,checkIntegerish)
-importFrom(checkmate,checkNumber)
-importFrom(checkmate,checkNumeric)
-importFrom(data.table,.I)
-importFrom(data.table,between)
-importFrom(data.table,data.table)
-importFrom(ggplot2,aes_string)
-importFrom(ggplot2,coord_fixed)
-importFrom(ggplot2,element_blank)
-importFrom(ggplot2,facet_wrap)
-importFrom(ggplot2,geom_abline)
-importFrom(ggplot2,geom_boxplot)
-importFrom(ggplot2,geom_hline)
-importFrom(ggplot2,geom_point)
-importFrom(ggplot2,geom_smooth)
-importFrom(ggplot2,geom_tile)
-importFrom(ggplot2,geom_violin)
-importFrom(ggplot2,ggplot)
-importFrom(ggplot2,ggtitle)
-importFrom(ggplot2,scale_colour_manual)
-importFrom(ggplot2,scale_fill_distiller)
-importFrom(ggplot2,scale_fill_manual)
-importFrom(ggplot2,scale_x_log10)
-importFrom(ggplot2,scale_y_continuous)
-importFrom(ggplot2,scale_y_log10)
-importFrom(ggplot2,theme)
-importFrom(ggplot2,theme_minimal)
-importFrom(ggplot2,xlab)
-importFrom(ggplot2,ylab)
-importFrom(locfit,locfit)
-importFrom(methods,"slot<-")
-importFrom(methods,as)
-importFrom(methods,callNextMethod)
-importFrom(methods,is)
-importFrom(methods,new)
-importFrom(methods,show)
-importFrom(methods,slot)
-importFrom(methods,slotNames)
-importFrom(methods,validObject)
-importFrom(stats,aggregate)
-importFrom(stats,approxfun)
-importFrom(stats,cor)
-importFrom(stats,dbeta)
-importFrom(stats,density)
-importFrom(stats,dnbinom)
-importFrom(stats,ks.test)
-importFrom(stats,median)
-importFrom(stats,model.matrix)
-importFrom(stats,nls)
-importFrom(stats,pnorm)
-importFrom(stats,qgamma)
-importFrom(stats,qnorm)
-importFrom(stats,quantile)
-importFrom(stats,rbinom)
-importFrom(stats,rchisq)
-importFrom(stats,rgamma)
-importFrom(stats,rlnorm)
-importFrom(stats,rnbinom)
-importFrom(stats,rnorm)
-importFrom(stats,rpois)
-importFrom(stats,runif)
importFrom(stats,sd)
-importFrom(stats,shapiro.test)
-importFrom(utils,globalVariables)
-importFrom(utils,head)
-importFrom(utils,read.delim)
-importFrom(utils,write.table)
diff --git a/R/data.R b/R/data.R
index 6f6a7a1..e69de29 100644
--- a/R/data.R
+++ b/R/data.R
@@ -1,76 +0,0 @@
-#' Example GFF file
-#'
-#' An example GFF file containing human genes from chromosome 22 of the GRCh38
-#'
-#' @docType data
-#'
-#' @usage data(ex_gff)
-#'
-#' @format A data frame with 504 rows and 9 variables:
-#' \describe{
-#' \item{V1}{id}
-#' \item{V2}{source}
-#' \item{V3}{feature}
-#' \item{V4}{start}
-#' \item{V5}{end}
-#' \item{V6}{score}
-#' \item{V7}{direction}
-#' \item{V8}{frame}
-#' \item{V9}{attribute}
-#' ...
-#' }
-#' @source \url{ftp://ftp.ensembl.org/pub/release-98/gff3/homo_sapiens/Homo_sapiens.GRCh38.98.chromosome.22.gff3.gz
-#'}
-"ex_gff"
-
-
-#' Example simulated genotype file (vcf)
-#'
-#' An example genotype file in the .vcf format, where each column is a sample
-#' and each row is a SNP.
-#'
-#' @docType data
-#'
-#' @usage data(ex_snps)
-#'
-#' @format A data frame with 141882 rows and 109 variables:
-#'
-#' @source \url{ftp://ftp.ensembl.org/pub/release-98/gff3/homo_sapiens/Homo_sapiens.GRCh38.98.chromosome.22.gff3.gz
-#'}
-"ex_snps"
-
-
-#' Example gene mean data for a whole population (from GTEx)
-#'
-#' An example population wide gene mean data.frame, where each column is a
-#' sample and each row is a gene. Data is from the Thyroid tissue from the GTEx
-#' dataset. Example data has already been filtered to remove genes with >50%
-#' of samples have < 0.1 TPM
-#'
-#' @docType data
-#'
-#' @usage data(ex_means)
-#'
-#' @format A data frame with 2438 rows and 850 variables:
-#'
-#' @source \url{https://storage.googleapis.com/gtex_analysis_v8/rna_seq_data/GTEx_Analysis_2017-06-05_v8_RNASeQCv1.1.9_gene_tpm.gct.gz
-#'}
-"ex_means"
-
-
-#' Example eQTL results from GTEx Thyroid tissue
-#'
-#' An example dataframe with eQTL pairs and effect sizes from the Thyroid
-#' tissue from the GTEx dataset. This example data has already been filtered to
-#' include only the top eSNP for each gene. Columns that must be included are:
-#' 'gene_id', 'pval_nominal', and 'slope'.
-#'
-#' @docType data
-#'
-#' @usage data(ex_pairs)
-#'
-#' @format A data frame with 25,244 rows and 9 variables:
-#'
-#' @source \url{https://storage.googleapis.com/gtex_analysis_v7/single_tissue_eqtl_data/all_snp_gene_associations/Thyroid.allpairs.txt.gz
-#'}
-"ex_pairs"
diff --git a/R/epi-functions.R b/R/epi-functions.R
new file mode 100644
index 0000000..2d37367
--- /dev/null
+++ b/R/epi-functions.R
@@ -0,0 +1,1524 @@
+#' Set control params
+#'
+#' Sets the controls for stochastic individual contact models simulated with
+#' icm.
+#'
+#' @param 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.
+#' @param 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.
+#' @param e.num Initial number of E compartment individuals in
+#' the simulated population.
+#' @param i.num Initial number of I compartment individuals in
+#' the simulated population.
+#' @param q.num Initial number of Q compartment individuals in
+#' the simulated population.
+#' @param h.num Initial number of H compartment individuals in
+#' the simulated population.
+#' @param r.num Initial number of R compartment individuals in
+#' the simulated population.
+#' @param f.num Initial number of F compartment individuals in
+#' the simulated population.
+#'
+#' @return control.icm object
+#'
+set.control <- function(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) {
+
+ control <- control.icm(type = type,
+ nsteps = nsteps,
+ nsims = nsims,
+ ncores = ncores,
+ prog.rand = prog.rand,
+ rec.rand = rec.rand,
+ fat.rand = fat.rand,
+ quar.rand = quar.rand,
+ hosp.rand = hosp.rand,
+ disch.rand = disch.rand)
+ return(control)
+}
+
+
+#' Set init params
+#'
+#' Sets the initial conditions for stochastic individual contact models
+#' simulated with icm.
+#'
+#' @param 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.
+#' @param e.num Initial number of E compartment individuals in
+#' the simulated population.
+#' @param i.num Initial number of I compartment individuals in
+#' the simulated population.
+#' @param q.num Initial number of Q compartment individuals in
+#' the simulated population.
+#' @param h.num Initial number of H compartment individuals in
+#' the simulated population.
+#' @param r.num Initial number of R compartment individuals in
+#' the simulated population.
+#' @param f.num Initial number of F compartment individuals in
+#' the simulated population.
+#'
+#' @return init.icm object
+#'
+set.init <- function(s.num = 9997,
+ e.num=0,
+ i.num = 3,
+ q.num=0,
+ h.num=0,
+ r.num = 0,
+ f.num = 0){
+
+ init <- init.icm(s.num = s.num,
+ e.num = e.num,
+ i.num = i.num,
+ q.num = q.num,
+ h.num = h.num,
+ r.num = r.num,
+ f.num = f.num)
+ return(init)
+}
+
+#' Set epidemic params
+#'
+#' Sets the epidemic parameters for stochastic individual contact models
+#' simulated with icm.
+#'
+#' @param 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.
+#' @param 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.
+#' @param 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.
+#' @param 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.
+#' @param 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.
+#' @param 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.
+#' @param 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.
+#' @param 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).
+#' @param disch.rate Rate per day at which people needing
+#' hospitalisation recover.
+#' @param 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.
+#' @param prog.dist.scale Scale parameter for Weibull distribution
+#' for progression, see prog.rand for details.
+#' @param 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.
+#' @param 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.
+#' @param rec.dist.scale Scale parameter for Weibull distribution for
+#' recovery, see rec.rand for details.
+#' @param 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.
+#' @param fat.rate.base Baseline mortality rate per day for people
+#' needing hospitalisation (deaths due to the virus). See fat.rand for more
+#' details.
+#' @param hosp.cap Number of available hospital beds for the modelled
+#' population. See fat.rand for more details.
+#' @param 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.
+#' @param 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.
+#' @param vital Enables demographics, that is, arrivals and
+#' departures, to and from the simulated population.
+#' @param 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.
+#' @param 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.
+#' @param 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.
+#' @param 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.
+#' @param 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.
+#' @param 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.
+#' @param 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.
+#' @param out Summary function for the simulation runs. median is
+#' also available, or percentiles, see the EpiModel documentation.
+#'
+#' @return param.icm object
+#'
+set.param <- function(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 = 1/20,
+ 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"){
+
+ param <- param.icm(inf.prob.e = inf.prob.e,
+ act.rate.e = act.rate.e,
+ inf.prob.i = inf.prob.i,
+ act.rate.i = act.rate.i,
+ inf.prob.q = inf.prob.q,
+ act.rate.q = act.rate.q,
+ quar.rate = quar.rate,
+ hosp.rate = hosp.rate,
+ disch.rate = disch.rate,
+ prog.rate = prog.rate,
+ prog.dist.scale = prog.dist.scale,
+ prog.dist.shape = prog.dist.shape,
+ rec.rate = rec.rate,
+ rec.dist.scale = rec.dist.scale,
+ rec.dist.shape = rec.dist.shape,
+ fat.rate.base = fat.rate.base,
+ hosp.cap = hosp.cap,
+ fat.rate.overcap = fat.rate.overcap,
+ fat.tcoeff = fat.tcoeff,
+ vital = vital,
+ a.rate = a.rate,
+ a.prop.e = a.prop.e,
+ a.prop.i = a.prop.i,
+ a.prop.q = a.prop.q,
+ ds.rate = ds.rate,
+ de.rate = de.rate,
+ di.rate = di.rate,
+ dq.rate = dq.rate,
+ dh.rate = dh.rate,
+ dr.rate = dr.rate)
+ return(param)
+}
+
+
+#' Initialize ICM
+#'
+#' Function to initialize the icm
+#'
+#' @param param ICM parameters.
+#' @param init Initial value parameters.
+#' @param control Control parameters.
+#'
+#' @return Updated dat
+initialize.FUN <- function(param, init, control) {
+
+ ## Master List for Data ##
+ dat <- list()
+ dat$param <- getParams(param, c(slotNames(param)))
+ dat$init <- getParams(init, c(slotNames(init)))
+ dat$control <- getParams(control, c(slotNames(control)))
+
+ # Set attributes
+ dat$attr <- list()
+ numeric.init <- init[which(sapply(init, class) == "numeric")]
+ n <- do.call("sum", numeric.init)
+ dat$attr$active <- rep(1, n)
+ if (dat$param$groups == 1) {
+ dat$attr$group <- rep(1, n)
+ } else {
+ g2inits <- grep(".g2", names(numeric.init))
+ g1inits <- setdiff(1:length(numeric.init), g2inits)
+ nG1 <- sum(sapply(g1inits, function(x) init[[x]]))
+ nG2 <- sum(sapply(g2inits, function(x) init[[x]]))
+ dat$attr$group <- c(rep(1, nG1), rep(2, max(0, nG2)))
+ }
+
+ # Initialize status and infection time
+ dat <- init_status.icm(dat)
+
+
+ # Summary out list
+ dat <- get_prev.icm(dat, at = 1)
+
+ return(dat)
+}
+
+
+#' Infection function
+#'
+#' Function to implement infection processes for SEIQHRF model
+#'
+#' @param dat Input data needed.
+#' @param at Not sure???
+#'
+#' @return Updated dat
+
+infection.FUN <- function(dat, at) {
+ type <- dat$control$type
+
+ # the following checks need to be moved to control.icm in due course
+ nsteps <- dat$control$nsteps
+
+ act.rate.i <- dat$param$act.rate.i
+ if (!(length(act.rate.i) == 1 || length(act.rate.i == nsteps))) {
+ stop("Length of act.rate.i must be 1 or the value of nsteps")
+ }
+ act.rate.i.g2 <- dat$param$act.rate.i.g2
+ if (!is.null(act.rate.i.g2) &&
+ !(length(act.rate.i.g2) == 1 || length(act.rate.i.g2 == nsteps))) {
+ stop("Length of act.rate.i.g2 must be 1 or the value of nsteps")
+ }
+ inf.prob.i <- dat$param$inf.prob.i
+ if (!(length(inf.prob.i) == 1 || length(inf.prob.i == nsteps))) {
+ stop("Length of inf.prob.i must be 1 or the value of nsteps")
+ }
+ inf.prob.i.g2 <- dat$param$inf.prob.i.g2
+ if (!is.null(inf.prob.i.g2) &&
+ !(length(inf.prob.i.g2) == 1 || length(inf.prob.i.g2 == nsteps))) {
+ stop("Length of inf.prob.i.g2 must be 1 or the value of nsteps")
+ }
+ if (type %in% c("SEIQHR", "SEIQHRF")) {
+ quar.rate <- dat$param$quar.rate
+ if (!(length(quar.rate) == 1 || length(quar.rate == nsteps))) {
+ stop("Length of quar.rate must be 1 or the value of nsteps")
+ }
+ quar.rate.g2 <- dat$param$quar.rate.g2
+ if (!is.null(quar.rate.g2) &&
+ !(length(quar.rate.g2) == 1 || length(quar.rate.g2 == nsteps))) {
+ stop("Length of quar.rate.g2 must be 1 or the value of nsteps")
+ }
+ disch.rate <- dat$param$disch.rate
+ if (!(length(disch.rate) == 1 || length(disch.rate == nsteps))) {
+ stop("Length of disch.rate must be 1 or the value of nsteps")
+ }
+ disch.rate.g2 <- dat$param$disch.rate.g2
+ if (!is.null(disch.rate.g2) &&
+ !(length(disch.rate.g2) == 1 || length(disch.rate.g2 == nsteps))) {
+ stop("Length of disch.rate.g2 must be 1 or the value of nsteps")
+ }
+ }
+
+ if (type %in% c("SEIQHR", "SEIQHRF")) {
+ act.rate.e <- dat$param$act.rate.e
+ if (!(length(act.rate.e) == 1 || length(act.rate.e == nsteps))) {
+ stop("Length of act.rate.e must be 1 or the value of nsteps")
+ }
+ act.rate.e.g2 <- dat$param$act.rate.e.g2
+ if (!is.null(act.rate.e.g2) &&
+ !(length(act.rate.e.g2) == 1 || length(act.rate.e.g2 == nsteps))) {
+ stop("Length of act.rate.e.g2 must be 1 or the value of nsteps")
+ }
+ inf.prob.e <- dat$param$inf.prob.e
+ if (!(length(inf.prob.e) == 1 || length(inf.prob.e == nsteps))) {
+ stop("Length of inf.prob.e must be 1 or the value of nsteps")
+ }
+ inf.prob.e.g2 <- dat$param$inf.prob.e.g2
+ if (!is.null(inf.prob.e.g2) &&
+ !(length(inf.prob.e.g2) == 1 || length(inf.prob.e.g2 == nsteps))) {
+ stop("Length of inf.prob.e.g2 must be 1 or the value of nsteps")
+ }
+
+ act.rate.q <- dat$param$act.rate.q
+ if (!(length(act.rate.q) == 1 || length(act.rate.q == nsteps))) {
+ stop("Length of act.rate.q must be 1 or the value of nsteps")
+ }
+ act.rate.q.g2 <- dat$param$act.rate.q.g2
+ if (!is.null(act.rate.q.g2) &&
+ !(length(act.rate.q.g2) == 1 || length(act.rate.q.g2 == nsteps))) {
+ stop("Length of act.rate.q.g2 must be 1 or the value of nsteps")
+ }
+ inf.prob.q <- dat$param$inf.prob.q
+ if (!(length(inf.prob.q) == 1 || length(inf.prob.q == nsteps))) {
+ stop("Length of inf.prob.q must be 1 or the value of nsteps")
+ }
+ inf.prob.q.g2 <- dat$param$inf.prob.q.g2
+ if (!is.null(inf.prob.q.g2) &&
+ !(length(inf.prob.q.g2) == 1 || length(inf.prob.q.g2 == nsteps))) {
+ stop("Length of inf.prob.q.g2 must be 1 or the value of nsteps")
+ }
+ }
+
+ # Transmission from infected
+ ## Expected acts
+ if (dat$param$groups == 1) {
+ if (length(act.rate.i) > 1) {
+ acts <- round(act.rate.i[at - 1] * dat$epi$num[at - 1] / 2)
+ } else {
+ acts <- round(act.rate.i * dat$epi$num[at - 1] / 2)
+ }
+ }
+ if (dat$param$groups == 2) {
+ if (dat$param$balance == "g1") {
+ if (length(act.rate.i) > 1) {
+ acts <- round(act.rate.i[at - 1] *
+ (dat$epi$num[at - 1] + dat$epi$num.g2[at - 1]) / 2)
+ } else {
+ acts <- round(act.rate.i *
+ (dat$epi$num[at - 1] + dat$epi$num.g2[at - 1]) / 2)
+ }
+ }
+ if (dat$param$balance == "g2") {
+ if (length(act.rate.i.g2) > 1) {
+ acts <- round(act.rate.i.g2[at - 1] *
+ (dat$epi$num[at - 1] + dat$epi$num.g2[at - 1]) / 2)
+ } else {
+ acts <- round(act.rate.i.g2 *
+ (dat$epi$num[at - 1] + dat$epi$num.g2[at - 1]) / 2)
+ }
+ }
+ }
+
+ ## Edgelist
+ if (dat$param$groups == 1) {
+ p1 <- ssample(which(dat$attr$active == 1 & dat$attr$status != "f"), acts, replace = TRUE)
+ p2 <- ssample(which(dat$attr$active == 1 & dat$attr$status != "f"), acts, replace = TRUE)
+ } else {
+ p1 <- ssample(which(dat$attr$active == 1 & dat$attr$group == 1 & dat$attr$status != "f"),
+ acts, replace = TRUE)
+ p2 <- ssample(which(dat$attr$active == 1 & dat$attr$group == 2 & dat$attr$status != "f"),
+ acts, replace = TRUE)
+ }
+
+ del <- NULL
+ if (length(p1) > 0 & length(p2) > 0) {
+ del <- data.frame(p1, p2)
+ if (dat$param$groups == 1) {
+ while (any(del$p1 == del$p2)) {
+ del$p2 <- ifelse(del$p1 == del$p2,
+ ssample(which(dat$attr$active == 1 & dat$attr$status != "f"), 1), del$p2)
+ }
+ }
+ }
+
+ ## Discordant edgelist (del)
+ del$p1.stat <- dat$attr$status[del$p1]
+ del$p2.stat <- dat$attr$status[del$p2]
+ serodis <- (del$p1.stat == "s" & del$p2.stat == "i") |
+ (del$p1.stat == "i" & del$p2.stat == "s")
+ del <- del[serodis == TRUE, ]
+
+ ## Transmission on edgelist
+ if (nrow(del) > 0) {
+ if (dat$param$groups == 1) {
+ if (length(inf.prob.i) > 1) {
+ del$tprob <- inf.prob.i[at]
+ } else {
+ del$tprob <- inf.prob.i
+ }
+ } else {
+ if (length(inf.prob.i) > 1) {
+ del$tprob <- ifelse(del$p1.stat == "s", inf.prob.i[at],
+ inf.prob.i.g2[at])
+ } else {
+ del$tprob <- ifelse(del$p1.stat == "s", inf.prob.i,
+ inf.prob.i.g2)
+ }
+ }
+ if (!is.null(dat$param$inter.eff.i) && at >= dat$param$inter.start.i &&
+ at <= dat$param$inter.stop.i) {
+ del$tprob <- del$tprob * (1 - dat$param$inter.eff.i)
+ }
+ del$trans <- rbinom(nrow(del), 1, del$tprob)
+ del <- del[del$trans == TRUE, ]
+ if (nrow(del) > 0) {
+ if (dat$param$groups == 1) {
+ newIds <- unique(ifelse(del$p1.stat == "s", del$p1, del$p2))
+ nExp.i <- length(newIds)
+ }
+ if (dat$param$groups == 2) {
+ newIdsg1 <- unique(del$p1[del$p1.stat == "s"])
+ newIdsg2 <- unique(del$p2[del$p2.stat == "s"])
+ nExp.i <- length(newIdsg1)
+ nExpg2.i <- length(newIdsg2)
+ newIds <- c(newIdsg1, newIdsg2)
+ }
+ dat$attr$status[newIds] <- "e"
+ dat$attr$expTime[newIds] <- at
+ } else {
+ nExp.i <- nExpg2.i <- 0
+ }
+ } else {
+ nExp.i <- nExpg2.i <- 0
+ }
+
+ if (type %in% c("SEIQHRF")) {
+
+ # Transmission from exposed
+ ## Expected acts
+ if (dat$param$groups == 1) {
+ if (length(act.rate.e) > 1) {
+ acts <- round(act.rate.e[at - 1] * dat$epi$num[at - 1] / 2)
+ } else {
+ acts <- round(act.rate.e * dat$epi$num[at - 1] / 2)
+ }
+ }
+ if (dat$param$groups == 2) {
+ if (dat$param$balance == "g1") {
+ if (length(act.rate.e.g2) > 1) {
+ acts <- round(act.rate.e.g2[at - 1] *
+ (dat$epi$num[at - 1] + dat$epi$num.g2[at - 1]) / 2)
+ } else {
+ acts <- round(act.rate.e.g2 *
+ (dat$epi$num[at - 1] + dat$epi$num.g2[at - 1]) / 2)
+ }
+ }
+ if (dat$param$balance == "g2") {
+ if (length(act.rate.e.g2) > 1) {
+ acts <- round(act.rate.e.g2[at - 1] *
+ (dat$epi$num[at - 1] + dat$epi$num.g2[at - 1]) / 2)
+ } else {
+ acts <- round(act.rate.e.g2 *
+ (dat$epi$num[at - 1] + dat$epi$num.g2[at - 1]) / 2)
+ }
+ }
+ }
+
+ ## Edgelist
+ if (dat$param$groups == 1) {
+ p1 <- ssample(which(dat$attr$active == 1 & dat$attr$status != "f"), acts, replace = TRUE)
+ p2 <- ssample(which(dat$attr$active == 1 & dat$attr$status != "f"), acts, replace = TRUE)
+ } else {
+ p1 <- ssample(which(dat$attr$active == 1 & dat$attr$group == 1 & dat$attr$status != "f"),
+ acts, replace = TRUE)
+ p2 <- ssample(which(dat$attr$active == 1 & dat$attr$group == 2 & dat$attr$status != "f"),
+ acts, replace = TRUE)
+ }
+
+ del <- NULL
+ if (length(p1) > 0 & length(p2) > 0) {
+ del <- data.frame(p1, p2)
+ if (dat$param$groups == 1) {
+ while (any(del$p1 == del$p2)) {
+ del$p2 <- ifelse(del$p1 == del$p2,
+ ssample(which(dat$attr$active == 1 & dat$attr$status != "f"), 1), del$p2)
+ }
+ }
+
+ ## Discordant edgelist (del)
+ del$p1.stat <- dat$attr$status[del$p1]
+ del$p2.stat <- dat$attr$status[del$p2]
+ # serodiscordance
+ serodis <- (del$p1.stat == "s" & del$p2.stat == "e") |
+ (del$p1.stat == "e" & del$p2.stat == "s")
+ del <- del[serodis == TRUE, ]
+
+ ## Transmission on edgelist
+ if (nrow(del) > 0) {
+ if (dat$param$groups == 1) {
+ if (length(inf.prob.e) > 1) {
+ del$tprob <- inf.prob.e[at]
+ } else {
+ del$tprob <- inf.prob.e
+ }
+ } else {
+ if (length(inf.prob.e) > 1) {
+ del$tprob <- ifelse(del$p1.stat == "s", inf.prob.e[at],
+ inf.prob.e.g2[at])
+ } else {
+ del$tprob <- ifelse(del$p1.stat == "s", inf.prob.e,
+ inf.prob.e.g2)
+ }
+ }
+ if (!is.null(dat$param$inter.eff.e) && at >= dat$param$inter.start.e &&
+ at <= dat$param$inter.stop.e) {
+ del$tprob <- del$tprob * (1 - dat$param$inter.eff.e)
+ }
+ del$trans <- rbinom(nrow(del), 1, del$tprob)
+ del <- del[del$trans == TRUE, ]
+ if (nrow(del) > 0) {
+ if (dat$param$groups == 1) {
+ newIds <- unique(ifelse(del$p1.stat == "s", del$p1, del$p2))
+ nExp.e <- length(newIds)
+ }
+ if (dat$param$groups == 2) {
+ newIdsg1 <- unique(del$p1[del$p1.stat == "s"])
+ newIdsg2 <- unique(del$p2[del$p2.stat == "s"])
+ nExp.e <- length(newIdsg1)
+ nExpg2.e <- length(newIdsg2)
+ newIds <- c(newIdsg1, newIdsg2)
+ }
+ dat$attr$status[newIds] <- "e"
+ dat$attr$expTime[newIds] <- at
+ } else {
+ nExp.e <- nExpg2.e <- 0
+ }
+ } else {
+ nExp.e <- nExpg2.e <- 0
+ }
+ } else {
+ nExp.e <- nExpg2.e <- 0
+ }
+
+
+ # Transmission from quarantined
+ ## Expected acts
+ if (dat$param$groups == 1) {
+ if (length(act.rate.q) > 1) {
+ acts <- round(act.rate.q[at - 1] * dat$epi$num[at - 1] / 2)
+ } else {
+ acts <- round(act.rate.q * dat$epi$num[at - 1] / 2)
+ }
+ }
+ if (dat$param$groups == 2) {
+ if (dat$param$balance == "g1") {
+ if (length(act.rate.q.g2) > 1) {
+ acts <- round(act.rate.q.g2[at - 1] *
+ (dat$epi$num[at - 1] + dat$epi$num.g2[at - 1]) / 2)
+ } else {
+ acts <- round(act.rate.q.g2 *
+ (dat$epi$num[at - 1] + dat$epi$num.g2[at - 1]) / 2)
+ }
+ }
+ if (dat$param$balance == "g2") {
+ if (length(act.rate.q.g2) > 1) {
+ acts <- round(act.rate.q.g2[at - 1] *
+ (dat$epi$num[at - 1] + dat$epi$num.g2[at - 1]) / 2)
+ } else {
+ acts <- round(act.rate.q.g2 *
+ (dat$epi$num[at - 1] + dat$epi$num.g2[at - 1]) / 2)
+ }
+ }
+ }
+
+ ## Edgelist
+ if (dat$param$groups == 1) {
+ p1 <- ssample(which(dat$attr$active == 1 & dat$attr$status != "f"), acts, replace = TRUE)
+ p2 <- ssample(which(dat$attr$active == 1 & dat$attr$status != "f"), acts, replace = TRUE)
+ } else {
+ p1 <- ssample(which(dat$attr$active == 1 & dat$attr$group == 1 & dat$attr$status != "f"),
+ acts, replace = TRUE)
+ p2 <- ssample(which(dat$attr$active == 1 & dat$attr$group == 2 & dat$attr$status != "f"),
+ acts, replace = TRUE)
+ }
+
+ del <- NULL
+ if (length(p1) > 0 & length(p2) > 0) {
+ del <- data.frame(p1, p2)
+ if (dat$param$groups == 1) {
+ while (any(del$p1 == del$p2)) {
+ del$p2 <- ifelse(del$p1 == del$p2,
+ ssample(which(dat$attr$active == 1 & dat$attr$status != "f"), 1), del$p2)
+ }
+ }
+
+ ## Discordant edgelist (del)
+ del$p1.stat <- dat$attr$status[del$p1]
+ del$p2.stat <- dat$attr$status[del$p2]
+ # serodiscordance
+ serodis <- (del$p1.stat == "s" & del$p2.stat == "q") |
+ (del$p1.stat == "q" & del$p2.stat == "s")
+ del <- del[serodis == TRUE, ]
+
+ ## Transmission on edgelist
+ if (nrow(del) > 0) {
+ if (dat$param$groups == 1) {
+ if (length(inf.prob.q) > 1) {
+ del$tprob <- inf.prob.q[at]
+ } else {
+ del$tprob <- inf.prob.q
+ }
+ } else {
+ if (length(inf.prob.q) > 1) {
+ del$tprob <- ifelse(del$p1.stat == "s", inf.prob.q[at],
+ inf.prob.q.g2[at])
+ } else {
+ del$tprob <- ifelse(del$p1.stat == "s", inf.prob.q,
+ inf.prob.q.g2)
+ }
+ }
+ if (!is.null(dat$param$inter.eff.q) && at >= dat$param$inter.start.q &&
+ at <= dat$param$inter.stop.q) {
+ del$tprob <- del$tprob * (1 - dat$param$inter.eff.q)
+ }
+ del$trans <- rbinom(nrow(del), 1, del$tprob)
+ del <- del[del$trans == TRUE, ]
+ if (nrow(del) > 0) {
+ if (dat$param$groups == 1) {
+ newIds <- unique(ifelse(del$p1.stat == "s", del$p1, del$p2))
+ nExp.q <- length(newIds)
+ }
+ if (dat$param$groups == 2) {
+ newIdsg1 <- unique(del$p1[del$p1.stat == "s"])
+ newIdsg2 <- unique(del$p2[del$p2.stat == "s"])
+ nExp.q <- length(newIdsg1)
+ nExpg2.q <- length(newIdsg2)
+ newIds <- c(newIdsg1, newIdsg2)
+ }
+ dat$attr$status[newIds] <- "e"
+ dat$attr$expTime[newIds] <- at
+ } else {
+ nExp.q <- nExpg2.q <- 0
+ }
+ } else {
+ nExp.q <- nExpg2.q <- 0
+ }
+ } else {
+ nExp.q <- nExpg2.q <- 0
+ }
+ }
+
+ ## Output
+ if (type %in% c("SEIQHR", "SEIQHRF")) {
+ if (at == 2) {
+ dat$epi$se.flow <- c(0, nExp.i + nExp.q)
+ } else {
+ dat$epi$se.flow[at] <- nExp.i + nExp.q
+ }
+ if (dat$param$groups == 2) {
+ if (at == 2) {
+ dat$epi$se.flow.g2 <- c(0, nExpg2.i + nExpg2.q )
+ } else {
+ dat$epi$se.flow.g2[at] <- nExpg2.i + nExpg2.q
+ }
+ }
+ } else {
+ if (at == 2) {
+ dat$epi$se.flow <- c(0, nExp.i)
+ } else {
+ dat$epi$se.flow[at] <- nExp.i
+ }
+ if (dat$param$groups == 2) {
+ if (at == 2) {
+ dat$epi$se.flow.g2 <- c(0, nExpg2.i)
+ } else {
+ dat$epi$se.flow.g2[at] <- nExpg2.i
+ }
+ }
+ }
+ return(dat)
+
+}
+
+
+#' Departures function
+#'
+#' Function to handel background demographics for the SEIQHRF model.
+#' Specifically departures (deaths not due to the virus, and emigration).
+#'
+#' @param dat Input data needed.
+#' @param at Not sure???
+#'
+#' @return Updated dat
+#'
+departures.FUN <- function(dat, at) {
+
+ # Conditions --------------------------------------------------------------
+ if (dat$param$vital == FALSE) {
+ return(dat)
+ }
+
+
+ # Variables ---------------------------------------------------------------
+ groups <- dat$param$groups
+ group <- dat$attr$group
+
+
+ # Susceptible departures ------------------------------------------------------
+ nDepartures <- nDeparturesG2 <- 0
+ idsElig <- which(dat$attr$active == 1 & dat$attr$status == "s")
+ nElig <- length(idsElig)
+ if (nElig > 0) {
+
+ gElig <- group[idsElig]
+ rates <- c(dat$param$ds.rate, dat$param$ds.rate.g2)
+ ratesElig <- rates[gElig]
+
+ if (dat$control$d.rand == TRUE) {
+ vecDepartures <- which(rbinom(nElig, 1, ratesElig) == 1)
+ if (length(vecDepartures) > 0) {
+ idsDpt <- idsElig[vecDepartures]
+ nDepartures <- sum(group[idsDpt] == 1)
+ nDeparturesG2 <- sum(group[idsDpt] == 2)
+ dat$attr$active[idsDpt] <- 0
+ }
+ } else {
+ nDepartures <- min(round(sum(ratesElig[gElig == 1])), sum(gElig == 1))
+ dat$attr$active[ssample(idsElig[gElig == 1], nDepartures)] <- 0
+ if (groups == 2) {
+ nDeparturesG2 <- min(round(sum(ratesElig[gElig == 2])), sum(gElig == 2))
+ dat$attr$active[ssample(idsElig[gElig == 2], nDeparturesG2)] <- 0
+ }
+ }
+ }
+
+ if (at == 2) {
+ dat$epi$ds.flow <- c(0, nDepartures)
+ if (groups == 2) {
+ dat$epi$ds.flow.g2 <- c(0, nDeparturesG2)
+ }
+ } else {
+ dat$epi$ds.flow[at] <- nDepartures
+ if (groups == 2) {
+ dat$epi$ds.flow.g2[at] <- nDeparturesG2
+ }
+ }
+
+ # Exposed Departures ---------------------------------------------------------
+ nDepartures <- nDeparturesG2 <- 0
+ idsElig <- which(dat$attr$active == 1 & dat$attr$status == "e")
+ nElig <- length(idsElig)
+ if (nElig > 0) {
+
+ gElig <- group[idsElig]
+ rates <- c(dat$param$de.rate, dat$param$de.rate.g2)
+ ratesElig <- rates[gElig]
+
+ if (dat$control$d.rand == TRUE) {
+ vecDepartures <- which(rbinom(nElig, 1, ratesElig) == 1)
+ if (length(vecDepartures) > 0) {
+ idsDpt <- idsElig[vecDepartures]
+ nDepartures <- sum(group[idsDpt] == 1)
+ nDeparturesG2 <- sum(group[idsDpt] == 2)
+ dat$attr$active[idsDpt] <- 0
+ }
+ } else {
+ nDepartures <- min(round(sum(ratesElig[gElig == 1])), sum(gElig == 1))
+ dat$attr$active[ssample(idsElig[gElig == 1], nDepartures)] <- 0
+ if (groups == 2) {
+ nDeparturesG2 <- min(round(sum(ratesElig[gElig == 2])), sum(gElig == 2))
+ dat$attr$active[ssample(idsElig[gElig == 2], nDeparturesG2)] <- 0
+ }
+ }
+ }
+
+ if (at == 2) {
+ dat$epi$de.flow <- c(0, nDepartures)
+ if (groups == 2) {
+ dat$epi$de.flow.g2 <- c(0, nDeparturesG2)
+ }
+ } else {
+ dat$epi$de.flow[at] <- nDepartures
+ if (groups == 2) {
+ dat$epi$de.flow.g2[at] <- nDeparturesG2
+ }
+ }
+
+
+ # Infected Departures ---------------------------------------------------------
+ nDepartures <- nDeparturesG2 <- 0
+ idsElig <- which(dat$attr$active == 1 & dat$attr$status == "i")
+ nElig <- length(idsElig)
+ if (nElig > 0) {
+
+ gElig <- group[idsElig]
+ rates <- c(dat$param$di.rate, dat$param$di.rate.g2)
+ ratesElig <- rates[gElig]
+
+ if (dat$control$d.rand == TRUE) {
+ vecDepartures <- which(rbinom(nElig, 1, ratesElig) == 1)
+ if (length(vecDepartures) > 0) {
+ idsDpt <- idsElig[vecDepartures]
+ nDepartures <- sum(group[idsDpt] == 1)
+ nDeparturesG2 <- sum(group[idsDpt] == 2)
+ dat$attr$active[idsDpt] <- 0
+ }
+ } else {
+ nDepartures <- min(round(sum(ratesElig[gElig == 1])), sum(gElig == 1))
+ dat$attr$active[ssample(idsElig[gElig == 1], nDepartures)] <- 0
+ if (groups == 2) {
+ nDeparturesG2 <- min(round(sum(ratesElig[gElig == 2])), sum(gElig == 2))
+ dat$attr$active[ssample(idsElig[gElig == 2], nDeparturesG2)] <- 0
+ }
+ }
+ }
+
+ if (at == 2) {
+ dat$epi$di.flow <- c(0, nDepartures)
+ if (groups == 2) {
+ dat$epi$di.flow.g2 <- c(0, nDeparturesG2)
+ }
+ } else {
+ dat$epi$di.flow[at] <- nDepartures
+ if (groups == 2) {
+ dat$epi$di.flow.g2[at] <- nDeparturesG2
+ }
+ }
+
+ # Quarantined Departures ---------------------------------------------------------
+ nDepartures <- nDeparturesG2 <- 0
+ idsElig <- which(dat$attr$active == 1 & dat$attr$status == "q")
+ nElig <- length(idsElig)
+ if (nElig > 0) {
+
+ gElig <- group[idsElig]
+ rates <- c(dat$param$dq.rate, dat$param$dq.rate.g2)
+ ratesElig <- rates[gElig]
+
+ if (dat$control$d.rand == TRUE) {
+ vecDepartures <- which(rbinom(nElig, 1, ratesElig) == 1)
+ if (length(vecDepartures) > 0) {
+ idsDpt <- idsElig[vecDepartures]
+ nDepartures <- sum(group[idsDpt] == 1)
+ nDeparturesG2 <- sum(group[idsDpt] == 2)
+ dat$attr$active[idsDpt] <- 0
+ }
+ } else {
+ nDepartures <- min(round(sum(ratesElig[gElig == 1])), sum(gElig == 1))
+ dat$attr$active[ssample(idsElig[gElig == 1], nDepartures)] <- 0
+ if (groups == 2) {
+ nDeparturesG2 <- min(round(sum(ratesElig[gElig == 2])), sum(gElig == 2))
+ dat$attr$active[ssample(idsElig[gElig == 2], nDeparturesG2)] <- 0
+ }
+ }
+ }
+
+ if (at == 2) {
+ dat$epi$dq.flow <- c(0, nDepartures)
+ if (groups == 2) {
+ dat$epi$dq.flow.g2 <- c(0, nDeparturesG2)
+ }
+ } else {
+ dat$epi$dq.flow[at] <- nDepartures
+ if (groups == 2) {
+ dat$epi$dq.flow.g2[at] <- nDeparturesG2
+ }
+ }
+
+ # Hospitalised Departures ---------------------------------------------------------
+ nDepartures <- nDeparturesG2 <- 0
+ idsElig <- which(dat$attr$active == 1 & dat$attr$status == "h")
+ nElig <- length(idsElig)
+ if (nElig > 0) {
+
+ gElig <- group[idsElig]
+ rates <- c(dat$param$dh.rate, dat$param$dh.rate.g2)
+ ratesElig <- rates[gElig]
+
+ if (dat$control$d.rand == TRUE) {
+ vecDepartures <- which(rbinom(nElig, 1, ratesElig) == 1)
+ if (length(vecDepartures) > 0) {
+ idsDpt <- idsElig[vecDepartures]
+ nDepartures <- sum(group[idsDpt] == 1)
+ nDeparturesG2 <- sum(group[idsDpt] == 2)
+ dat$attr$active[idsDpt] <- 0
+ }
+ } else {
+ nDepartures <- min(round(sum(ratesElig[gElig == 1])), sum(gElig == 1))
+ dat$attr$active[ssample(idsElig[gElig == 1], nDepartures)] <- 0
+ if (groups == 2) {
+ nDeparturesG2 <- min(round(sum(ratesElig[gElig == 2])), sum(gElig == 2))
+ dat$attr$active[ssample(idsElig[gElig == 2], nDeparturesG2)] <- 0
+ }
+ }
+ }
+
+ if (at == 2) {
+ dat$epi$dh.flow <- c(0, nDepartures)
+ if (groups == 2) {
+ dat$epi$dh.flow.g2 <- c(0, nDeparturesG2)
+ }
+ } else {
+ dat$epi$dh.flow[at] <- nDepartures
+ if (groups == 2) {
+ dat$epi$dh.flow.g2[at] <- nDeparturesG2
+ }
+ }
+
+
+ # Recovered Departures --------------------------------------------------------
+ nDepartures <- nDeparturesG2 <- 0
+ idsElig <- which(dat$attr$active == 1 & dat$attr$status == "r")
+ nElig <- length(idsElig)
+ if (nElig > 0) {
+
+ gElig <- group[idsElig]
+ rates <- c(dat$param$dr.rate, dat$param$dr.rate.g2)
+ ratesElig <- rates[gElig]
+
+ if (dat$control$d.rand == TRUE) {
+ vecDepartures <- which(rbinom(nElig, 1, ratesElig) == 1)
+ if (length(vecDepartures) > 0) {
+ idsDpt <- idsElig[vecDepartures]
+ nDepartures <- sum(group[idsDpt] == 1)
+ nDeparturesG2 <- sum(group[idsDpt] == 2)
+ dat$attr$active[idsDpt] <- 0
+ }
+ } else {
+ nDepartures <- min(round(sum(ratesElig[gElig == 1])), sum(gElig == 1))
+ dat$attr$active[ssample(idsElig[gElig == 1], nDepartures)] <- 0
+ if (groups == 2) {
+ nDeparturesG2 <- min(round(sum(ratesElig[gElig == 2])), sum(gElig == 2))
+ dat$attr$active[ssample(idsElig[gElig == 2], nDeparturesG2)] <- 0
+ }
+ }
+ }
+
+ if (at == 2) {
+ dat$epi$dr.flow <- c(0, nDepartures)
+ if (groups == 2) {
+ dat$epi$dr.flow.g2 <- c(0, nDeparturesG2)
+ }
+ } else {
+ dat$epi$dr.flow[at] <- nDepartures
+ if (groups == 2) {
+ dat$epi$dr.flow.g2[at] <- nDeparturesG2
+ }
+ }
+
+ return(dat)
+}
+
+
+#' Arrivals function
+#'
+#' Function to handel background demographics for the SEIQHRF model.
+#' 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..
+#'
+#' @param dat Input data needed.
+#' @param at Not sure???
+#'
+#' @return Updated dat
+#'
+arrivals.FUN <- function(dat, at) {
+
+ # Conditions --------------------------------------------------------------
+ if (dat$param$vital == FALSE) {
+ return(dat)
+ }
+
+ # Variables ---------------------------------------------------------------
+ a.rate <- dat$param$a.rate
+ a.rate.g2 <- dat$param$a.rate.g2
+ a.rand <- dat$control$a.rand
+ groups <- dat$param$groups
+ nOld <- dat$epi$num[at - 1]
+
+ # checking params, should be in control.icm or params.icm eventually
+ type <- dat$control$type
+ nsteps <- dat$control$nsteps
+
+ if (!(length(a.rate) == 1 || length(a.rate == nsteps))) {
+ stop("Length of a.rate must be 1 or the value of nsteps")
+ }
+ if (!is.null(a.rate.g2) &&
+ !(length(a.rate.g2) == 1 || length(a.rate.g2 == nsteps))) {
+ stop("Length of a.rate.g2 must be 1 or the value of nsteps")
+ }
+
+ a.prop.e <- dat$param$a.prop.e
+ if (!(length(a.prop.e) == 1 || length(a.prop.e == nsteps))) {
+ stop("Length of a.prop.e must be 1 or the value of nsteps")
+ }
+ a.prop.i <- dat$param$a.prop.i
+ if (!(length(a.prop.i) == 1 || length(a.prop.i == nsteps))) {
+ stop("Length of a.prop.i must be 1 or the value of nsteps")
+ }
+ a.prop.q <- dat$param$a.prop.q
+ if (!(length(a.prop.q) == 1 || length(a.prop.q == nsteps))) {
+ stop("Length of a.prop.q must be 1 or the value of nsteps")
+ }
+
+ a.prop.e.g2 <- dat$param$a.prop.e.g2
+ if (!is.null(a.prop.e.g2) &&
+ !(length(a.prop.e.g2) == 1 || length(a.prop.e.g2 == nsteps))) {
+ stop("Length of a.prop.e.g2 must be 1 or the value of nsteps")
+ }
+ a.prop.i.g2 <- dat$param$a.prop.i.g2
+ if (!is.null(a.prop.i.g2) &&
+ !(length(a.prop.i.g2) == 1 || length(a.prop.i.g2 == nsteps))) {
+ stop("Length of a.prop.i.g2 must be 1 or the value of nsteps")
+ }
+ a.prop.q.g2 <- dat$param$a.prop.q.g2
+ if (!is.null(a.prop.q.g2) &&
+ !(length(a.prop.q.g2) == 1 || length(a.prop.q.g2 == nsteps))) {
+ stop("Length of a.prop.q.g2 must be 1 or the value of nsteps")
+ }
+
+ # Process -----------------------------------------------------------------
+ nArrivals <- nArrivals.g2 <- 0
+
+ if (groups == 1) {
+ if (a.rand == TRUE) {
+ nArrivals <- sum(rbinom(nOld, 1, a.rate))
+ }
+ if (a.rand == FALSE) {
+ nArrivals <- round(nOld * a.rate)
+ }
+ }
+ if (groups == 2) {
+ nOldG2 <- dat$epi$num.g2[at - 1]
+ if (a.rand == TRUE) {
+ if (is.na(a.rate.g2)) {
+ nArrivals <- sum(rbinom(nOld, 1, a.rate))
+ nArrivals.g2 <- sum(rbinom(nOld, 1, a.rate))
+ } else {
+ nArrivals <- sum(rbinom(nOld, 1, a.rate))
+ nArrivals.g2 <- sum(rbinom(nOldG2, 1, a.rate.g2))
+ }
+ }
+ if (a.rand == FALSE) {
+ if (is.na(a.rate.g2)) {
+ nArrivals <- round(nOld * a.rate)
+ nArrivals.g2 <- round(nOld * a.rate)
+ } else {
+ nArrivals <- round(nOld * a.rate)
+ nArrivals.g2 <- round(nOldG2 * a.rate.g2)
+ }
+ }
+ }
+
+
+ ## Set attributes
+ totArrivals <- 0
+ totArrivals.g2 <- 0
+
+ # partition arrivals into compartments
+ if (length(a.prop.e) > 1) {
+ nArrivals.e <- round(nArrivals*(a.prop.e[at]))
+ totArrivals <- totArrivals + nArrivals.e
+ if (!is.null(a.prop.e.g2)) {
+ nArrivals.e.g2 <- round(nArrivals.g2*(a.prop.e.g2[at]))
+ totArrivals.g2 <- totArrivals.g2 + nArrivals.e.g2
+ } else {
+ nArrivals.e.g2 <- round(nArrivals.g2*(a.prop.e.g2[at]))
+ totArrivals.g2 <- totArrivals.g2 + nArrivals.e.g2
+ }
+ } else {
+ nArrivals.e <- round(nArrivals*a.prop.e)
+ totArrivals <- totArrivals + nArrivals.e
+ if (!is.null(a.prop.e.g2)) {
+ nArrivals.e.g2 <- round(nArrivals.g2*(a.prop.e.g2))
+ totArrivals.g2 <- totArrivals.g2 + nArrivals.e.g2
+ } else {
+ nArrivals.e.g2 <- round(nArrivals.g2*(a.prop.e.g2))
+ totArrivals.g2 <- totArrivals.g2 + nArrivals.e.g2
+ }
+ }
+
+ if (length(a.prop.i) > 1) {
+ nArrivals.i <- round(nArrivals*(a.prop.i[at]))
+ totArrivals <- totArrivals + nArrivals.i
+ if (!is.null(a.prop.i.g2)) {
+ nArrivals.i.g2 <- round(nArrivals.g2*(a.prop.i.g2[at]))
+ totArrivals.g2 <- totArrivals.g2 + nArrivals.i.g2
+ } else {
+ nArrivals.i.g2 <- round(nArrivals.g2*(a.prop.i.g2[at]))
+ totArrivals.g2 <- totArrivals.g2 + nArrivals.i.g2
+ }
+ } else {
+ nArrivals.i <- round(nArrivals*a.prop.i)
+ totArrivals <- totArrivals + nArrivals.i
+ if (!is.null(a.prop.i.g2)) {
+ nArrivals.i.g2 <- round(nArrivals.g2*(a.prop.i.g2))
+ totArrivals.g2 <- totArrivals.g2 + nArrivals.i.g2
+ } else {
+ nArrivals.i.g2 <- round(nArrivals.g2*(a.prop.i.g2))
+ totArrivals.g2 <- totArrivals.g2 + nArrivals.i.g2
+ }
+ }
+
+ if (length(a.prop.q) > 1) {
+ nArrivals.q <- round(nArrivals*(a.prop.q[at]))
+ totArrivals <- totArrivals + nArrivals.q
+ if (!is.null(a.prop.q.g2)) {
+ nArrivals.q.g2 <- round(nArrivals.g2*(a.prop.q.g2[at]))
+ totArrivals.g2 <- totArrivals.g2 + nArrivals.q.g2
+ } else {
+ nArrivals.q.g2 <- round(nArrivals.g2*(a.prop.q.g2[at]))
+ totArrivals.g2 <- totArrivals.g2 + nArrivals.q.g2
+ }
+ } else {
+ nArrivals.q <- round(nArrivals*a.prop.q)
+ totArrivals <- totArrivals + nArrivals.q
+ if (!is.null(a.prop.q.g2)) {
+ nArrivals.q.g2 <- round(nArrivals.g2*(a.prop.q.g2))
+ totArrivals.g2 <- totArrivals.g2 + nArrivals.q.g2
+ } else {
+ nArrivals.q.g2 <- round(nArrivals.g2*(a.prop.q.g2))
+ totArrivals.g2 <- totArrivals.g2 + nArrivals.q.g2
+ }
+ }
+
+ # debug
+ print("totArrivals:")
+ print(totArrivals)
+ print("totArrivals.g2:")
+ print(totArrivals.g2)
+ print("----")
+
+ # group 1
+ dat$attr$active <- c(dat$attr$active, rep(1, totArrivals))
+ dat$attr$group <- c(dat$attr$group, rep(1, totArrivals))
+ dat$attr$status <- c(dat$attr$status,
+ rep("e", nArrivals.e),
+ rep("i", nArrivals.i),
+ rep("q", nArrivals.q),
+ rep("s", totArrivals - nArrivals.e - nArrivals.i - nArrivals.q))
+ dat$attr$expTime <- c(dat$attr$expTime, rep(NA, totArrivals))
+ dat$attr$infTime <- c(dat$attr$infTime, rep(NA, totArrivals))
+ dat$attr$quarTime <- c(dat$attr$quarTime, rep(NA, totArrivals))
+ dat$attr$hospTime <- c(dat$attr$ihospTime, rep(NA, totArrivals))
+ dat$attr$recovTime <- c(dat$attr$recovTime, rep(NA, totArrivals))
+ dat$attr$fatTime <- c(dat$attr$fatTime, rep(NA, totArrivals))
+
+ # group 2
+ if (length(totArrivals.g2) > 0) {
+ dat$attr$active <- c(dat$attr$active, rep(1, totArrivals.g2))
+ dat$attr$group <- c(dat$attr$group, rep(2, totArrivals.g2))
+ dat$attr$status <- c(dat$attr$status,
+ rep("e", nArrivals.e.g2),
+ rep("i", nArrivals.i.g2),
+ rep("q", nArrivals.q.g2),
+ rep("s", totArrivals.g2 - nArrivals.e.g2 -
+ nArrivals.i.g2 - nArrivals.q.g2))
+ dat$attr$expTime <- c(dat$attr$expTime, rep(NA, totArrivals.g2))
+ dat$attr$infTime <- c(dat$attr$infTime, rep(NA, totArrivals.g2))
+ dat$attr$quarTime <- c(dat$attr$quarTime, rep(NA, totArrivals.g2))
+ dat$attr$hospTime <- c(dat$attr$ihospTime, rep(NA, totArrivals.g2))
+ dat$attr$recovTime <- c(dat$attr$recovTime, rep(NA, totArrivals.g2))
+ dat$attr$fatTime <- c(dat$attr$fatTime, rep(NA, totArrivals.g2))
+ }
+
+ # Output ------------------------------------------------------------------
+ if (at == 2) {
+ dat$epi$a.flow <- c(0, totArrivals)
+ dat$epi$a.e.flow <- c(0, nArrivals.e)
+ dat$epi$a.i.flow <- c(0, nArrivals.i)
+ dat$epi$a.q.flow <- c(0, nArrivals.q)
+ } else {
+ dat$epi$a.flow[at] <- totArrivals
+ dat$epi$a.e.flow[at] <- c(0, nArrivals.e)
+ dat$epi$a.i.flow[at] <- c(0, nArrivals.i)
+ dat$epi$a.q.flow[at] <- c(0, nArrivals.q)
+ }
+ if (length(totArrivals.g2) > 0 && dat$param$groups == 2) {
+ if (at == 2) {
+ dat$epi$a.flow.g2 <- c(0, totArrivals.g2)
+ dat$epi$a.e.flow.g2 <- c(0, nArrivals.e.g2)
+ dat$epi$a.i.flow.g2 <- c(0, nArrivals.i.g2)
+ dat$epi$a.q.flow.g2 <- c(0, nArrivals.q.g2)
+ } else {
+ dat$epi$a.flow.g2[at] <- totArrivals.g2
+ dat$epi$a.e.flow.g2[at] <- c(0, nArrivals.e.g2)
+ dat$epi$a.i.flow.g2[at] <- c(0, nArrivals.i.g2)
+ dat$epi$a.q.flow.g2[at] <- c(0, nArrivals.q.g2)
+ }
+ }
+
+ return(dat)
+}
+
+
+#' Get previous function
+#'
+#' Utility function that collects prevalence and transition time data from each
+#' run and stores it away in the simulation result object.
+#'
+#' @param dat Input data needed.
+#' @param at Not sure???
+#'
+#' @return Updated dat
+#'
+get_prev.FUN <- function(dat, at) {
+
+ if (at == 1) {
+ dat$epi <- list()
+ dat$epi$s.num <- sum(dat$attr$active == 1 &
+ dat$attr$status == "s" &
+ dat$attr$group == 1)
+ dat$epi$i.num <- sum(dat$attr$active == 1 &
+ dat$attr$status == "i" &
+ dat$attr$group == 1)
+ if (dat$control$type %in% c("SEIR", "SEIQHR", "SEIQHRF")) {
+ dat$epi$e.num <- sum(dat$attr$active == 1 &
+ dat$attr$status == "e" &
+ dat$attr$group == 1)
+ }
+ if (dat$control$type %in% c("SIR", "SEIR", "SEIQHR", "SEIQHRF")) {
+ dat$epi$r.num <- sum(dat$attr$active == 1 &
+ dat$attr$status == "r" &
+ dat$attr$group == 1)
+ }
+ if (dat$control$type %in% c("SEIQHR", "SEIQHRF")) {
+ dat$epi$q.num <- sum(dat$attr$active == 1 &
+ dat$attr$status == "q" &
+ dat$attr$group == 1)
+ dat$epi$h.num <- sum(dat$attr$active == 1 &
+ dat$attr$status == "h" &
+ dat$attr$group == 1)
+ }
+ if (dat$control$type =="SEIQHRF") {
+ dat$epi$f.num <- sum(dat$attr$active == 1 &
+ dat$attr$status == "f" &
+ dat$attr$group == 1)
+ }
+ if (dat$control$type == "SIR") {
+ dat$epi$num <- dat$epi$s.num +
+ dat$epi$i.num +
+ dat$epi$r.num
+ } else if (dat$control$type == "SEIR") {
+ dat$epi$num <- dat$epi$s.num +
+ dat$epi$e.num +
+ dat$epi$i.num +
+ dat$epi$r.num
+ } else if (dat$control$type == "SEIQHR") {
+ dat$epi$num <- dat$epi$s.num +
+ dat$epi$e.num +
+ dat$epi$i.num +
+ dat$epi$q.num +
+ dat$epi$h.num +
+ dat$epi$r.num
+ } else if (dat$control$type == "SEIQHRF") {
+ dat$epi$num <- dat$epi$s.num +
+ dat$epi$e.num +
+ dat$epi$i.num +
+ dat$epi$q.num +
+ dat$epi$h.num +
+ dat$epi$r.num +
+ dat$epi$f.num
+ } else {
+ dat$epi$num <- dat$epi$s.num + dat$epi$i.num
+ }
+ if (dat$param$groups == 2) {
+ dat$epi$s.num.g2 <- sum(dat$attr$active == 1 &
+ dat$attr$status == "s" &
+ dat$attr$group == 2)
+ if (dat$control$type %in% c("SEIR", "SEIQHR", "SEIQHRF")) {
+ dat$epi$e.num.g2 <- sum(dat$attr$active == 1 &
+ dat$attr$status == "e" &
+ dat$attr$group == 2)
+ }
+ if (dat$control$type %in% c("SEIQHR", "SEIQHRF")) {
+ dat$epi$q.num.g2 <- sum(dat$attr$active == 1 &
+ dat$attr$status == "q" &
+ dat$attr$group == 2)
+ dat$epi$h.num.g2 <- sum(dat$attr$active == 1 &
+ dat$attr$status == "h" &
+ dat$attr$group == 2)
+ }
+ if (dat$control$type %in% c("SEIQHRF")) {
+ dat$epi$f.num.g2 <- sum(dat$attr$active == 1 &
+ dat$attr$status == "f" &
+ dat$attr$group == 2)
+ }
+ dat$epi$i.num.g2 <- sum(dat$attr$active == 1 &
+ dat$attr$status == "i" &
+ dat$attr$group == 2)
+ dat$epi$num.g2 <- dat$epi$s.num.g2 + dat$epi$i.num.g2
+ if (dat$control$type %in% c("SIR", "SEIR", "SEIQHR", "SEIQHRF")) {
+ dat$epi$r.num.g2 <- sum(dat$attr$active == 1 &
+ dat$attr$status == "r" &
+ dat$attr$group == 2)
+ }
+ if (dat$control$type == "SIR") {
+ dat$epi$num.g2 <- dat$epi$s.num.g2 +
+ dat$epi$i.num.g2 +
+ dat$epi$r.num.g2
+ } else if (dat$control$type == "SEIR") {
+ dat$epi$num.g2 <- dat$epi$s.num.g2 +
+ dat$epi$e.num.g2 +
+ dat$epi$i.num.g2 +
+ dat$epi$r.num.g2
+ } else if (dat$control$type == "SEIQHR") {
+ dat$epi$num.g2 <- dat$epi$s.num.g2 +
+ dat$epi$e.num.g2 +
+ dat$epi$i.num.g2 +
+ dat$epi$q.num.g2 +
+ dat$epi$h.num.g2 +
+ dat$epi$r.num.g2
+ } else if (dat$control$type == "SEIQHRF") {
+ dat$epi$num.g2 <- dat$epi$s.num.g2 +
+ dat$epi$e.num.g2 +
+ dat$epi$i.num.g2 +
+ dat$epi$q.num.g2 +
+ dat$epi$h.num.g2 +
+ dat$epi$r.num.g2 +
+ dat$epi$f.num.g2
+ } else {
+ dat$epi$num.g2 <- dat$epi$s.num.g2 + dat$epi$i.num.g2
+ }
+ }
+ } else {
+ dat$epi$s.num[at] <- sum(dat$attr$active == 1 &
+ dat$attr$status == "s" &
+ dat$attr$group == 1)
+ if (dat$control$type %in% c("SEIR", "SEIQHR", "SEIQHRF")) {
+ dat$epi$e.num[at] <- sum(dat$attr$active == 1 &
+ dat$attr$status == "e" &
+ dat$attr$group == 1)
+ }
+ dat$epi$i.num[at] <- sum(dat$attr$active == 1 &
+ dat$attr$status == "i" &
+ dat$attr$group == 1)
+ if (dat$control$type %in% c("SIR", "SEIR", "SEIQHR", "SEIQHRF")) {
+ dat$epi$r.num[at] <- sum(dat$attr$active == 1 &
+ dat$attr$status == "r" &
+ dat$attr$group == 1)
+ }
+ if (dat$control$type %in% c("SEIQHR", "SEIQHRF")) {
+ dat$epi$q.num[at] <- sum(dat$attr$active == 1 &
+ dat$attr$status == "q" &
+ dat$attr$group == 1)
+ dat$epi$h.num[at] <- sum(dat$attr$active == 1 &
+ dat$attr$status == "h" &
+ dat$attr$group == 1)
+ }
+ if (dat$control$type %in% c("SEIQHRF")) {
+ dat$epi$f.num[at] <- sum(dat$attr$active == 1 &
+ dat$attr$status == "f" &
+ dat$attr$group == 1)
+ }
+ if (dat$control$type == "SIR") {
+ dat$epi$num[at] <- dat$epi$s.num[at] +
+ dat$epi$i.num[at] +
+ dat$epi$r.num[at]
+ } else if (dat$control$type == "SEIR") {
+ dat$epi$num[at] <- dat$epi$s.num[at] +
+ dat$epi$e.num[at] +
+ dat$epi$i.num[at] +
+ dat$epi$r.num[at]
+ } else if (dat$control$type == "SEIQHR") {
+ dat$epi$num[at] <- dat$epi$s.num[at] +
+ dat$epi$e.num[at] +
+ dat$epi$i.num[at] +
+ dat$epi$q.num[at] +
+ dat$epi$h.num[at] +
+ dat$epi$r.num[at]
+ } else if (dat$control$type == "SEIQHRF") {
+ dat$epi$num[at] <- dat$epi$s.num[at] +
+ dat$epi$e.num[at] +
+ dat$epi$i.num[at] +
+ dat$epi$q.num[at] +
+ dat$epi$h.num[at] +
+ dat$epi$r.num[at] +
+ dat$epi$f.num[at]
+ } else {
+ dat$epi$num[at] <- dat$epi$s.num[at] + dat$epi$i.num[at]
+ }
+
+ if (dat$param$groups == 2) {
+ dat$epi$s.num.g2[at] <- sum(dat$attr$active == 1 &
+ dat$attr$status == "s" &
+ dat$attr$group == 2)
+ if (dat$control$type %in% c("SEIR", "SEIQHR", "SEIQHRF")) {
+ dat$epi$i.num.g2[at] <- sum(dat$attr$active == 1 &
+ dat$attr$status == "e" &
+ dat$attr$group == 2)
+ }
+ dat$epi$i.num.g2[at] <- sum(dat$attr$active == 1 &
+ dat$attr$status == "i" &
+ dat$attr$group == 2)
+ if (dat$control$type %in% c("SIR", "SEIR", "SEIQHR", "SEIQHRF")) {
+ dat$epi$r.num.g2[at] <- sum(dat$attr$active == 1 &
+ dat$attr$status == "r" &
+ dat$attr$group == 2)
+ }
+ if (dat$control$type %in% c("SEIQHR", "SEIQHRF")) {
+ dat$epi$q.num.g2[at] <- sum(dat$attr$active == 1 &
+ dat$attr$status == "q" &
+ dat$attr$group == 2)
+ dat$epi$h.num.g2[at] <- sum(dat$attr$active == 1 &
+ dat$attr$status == "h" &
+ dat$attr$group == 2)
+ }
+ if (dat$control$type %in% c("SEIQHRF")) {
+ dat$epi$f.num.g2[at] <- sum(dat$attr$active == 1 &
+ dat$attr$status == "f" &
+ dat$attr$group == 2)
+ }
+ if (dat$control$type == "SIR") {
+ dat$epi$num.g2[at] <- dat$epi$s.num.g2[at] +
+ dat$epi$i.num.g2[at] +
+ dat$epi$r.num.g2[at]
+ } else if (dat$control$type == "SEIR") {
+ dat$epi$num.g2[at] <- dat$epi$s.num.g2[at] +
+ dat$epi$e.num.g2[at] +
+ dat$epi$i.num.g2[at] +
+ dat$epi$r.num.g2[at]
+ } else if (dat$control$type == "SEIQHR") {
+ dat$epi$num.g2[at] <- dat$epi$s.num.g2[at] +
+ dat$epi$e.num.g2[at] +
+ dat$epi$i.num.g2[at] +
+ dat$epi$q.num.g2[at] +
+ dat$epi$h.num.g2[at] +
+ dat$epi$r.num.g2[at]
+ } else if (dat$control$type == "SEIQHRF") {
+ dat$epi$num.g2[at] <- dat$epi$s.num.g2[at] +
+ dat$epi$e.num.g2[at] +
+ dat$epi$i.num.g2[at] +
+ dat$epi$q.num.g2[at] +
+ dat$epi$h.num.g2[at] +
+ dat$epi$r.num.g2[at] +
+ dat$epi$f.num.g2[at]
+ } else {
+ dat$epi$num.g2[at] <- dat$epi$s.num.g2[at] +
+ dat$epi$i.num.g2[at]
+ }
+ }
+ }
+
+ return(dat)
+}
+
+
diff --git a/R/results-parse.R b/R/results-parse.R
new file mode 100644
index 0000000..7c63959
--- /dev/null
+++ b/R/results-parse.R
@@ -0,0 +1,40 @@
+#' Get times
+#'
+#' Function to extract timings and assemble results in a dataframe
+#'
+#' @param simulate_results results from `simulte()` function.
+#'
+#' @return Dataframe of
+
+get_times <- function(simulate_results) {
+
+ sim <- simulate_results$sim
+
+ for (s in 1:sim$control$nsims) {
+ if (s == 1) {
+ times <- sim$times[[paste("sim", s, sep = "")]]
+ times <- times %>% mutate(s = s)
+ } else {
+ times <- times %>% bind_rows(sim$times[[paste("sim", s, sep = "")]]
+ %>% mutate(s = s))
+ }
+ }
+
+ times <- times %>% mutate(infTime = ifelse(infTime < 0, -5,
+ infTime), expTime = ifelse(expTime < 0, -5, expTime)) %>%
+ mutate(incubation_period = infTime - expTime, illness_duration = recovTime -
+ expTime, illness_duration_hosp = dischTime - expTime,
+ hosp_los = dischTime - hospTime, quarantine_delay = quarTime -
+ infTime, survival_time = fatTime - infTime) %>%
+ select(s, incubation_period, quarantine_delay, illness_duration,
+ illness_duration_hosp, hosp_los, survival_time) %>%
+ pivot_longer(-s, names_to = "period_type", values_to = "duration") %>%
+ mutate(period_type = factor(period_type, levels = c("incubation_period",
+ "quarantine_delay", "illness_duration", "illness_duration_hosp",
+ "hosp_los", "survival_time"), labels = c("Incubation period",
+ "Delay entering isolation", "Illness duration", "Illness duration (hosp)",
+ "Hospital care required duration", "Survival time of case fatalities"), ordered = TRUE))
+ return(times)
+}
+
+
diff --git a/R/sirPlus-functions.R b/R/sirPlus-functions.R
new file mode 100644
index 0000000..ec2dad1
--- /dev/null
+++ b/R/sirPlus-functions.R
@@ -0,0 +1,956 @@
+#' Infection function
+#'
+#' Function to implement infection processes for SEIQHRF model
+#'
+#' @param param ICM parameters.
+#' @param init Initial value parameters.
+#' @param control Control parameters.
+#'
+#' @return Updated dat
+#'
+icm.seiqhrf <- function(param, init, control) {
+
+ crosscheck.icm(param, init, control)
+ verbose.icm(control, type = "startup")
+ nsims <- control$nsims
+ ncores <- ifelse(control$nsims == 1, 1, min(future::availableCores(), control$ncores))
+ control$ncores <- ncores
+
+ if (ncores == 1) {
+
+ # Simulation loop start
+ for (s in 1:control$nsims) {
+
+ ## Initialization module
+ if (!is.null(control[["initialize.FUN"]])) {
+ dat <- do.call(control[["initialize.FUN"]], list(param, init, control))
+ }
+
+
+ # Timestep loop
+ for (at in 2:control$nsteps) {
+
+ ## User Modules
+ um <- control$user.mods
+ if (length(um) > 0) {
+ for (i in 1:length(um)) {
+ dat <- do.call(control[[um[i]]], list(dat, at))
+ }
+ }
+
+ ## Infection
+ if (!is.null(control[["infection.FUN"]])) {
+ dat <- do.call(control[["infection.FUN"]], list(dat, at))
+ }
+
+
+ ## Recovery
+ if (!is.null(control[["recovery.FUN"]])) {
+ dat <- do.call(control[["recovery.FUN"]], list(dat, at))
+ }
+
+
+ ## Departure Module
+ if (!is.null(control[["departures.FUN"]])) {
+ dat <- do.call(control[["departures.FUN"]], list(dat, at))
+ }
+
+
+ ## Arrival Module
+ if (!is.null(control[["arrivals.FUN"]])) {
+ dat <- do.call(control[["arrivals.FUN"]], list(dat, at))
+ }
+
+
+ ## Outputs
+ if (!is.null(control[["get_prev.FUN"]])) {
+ dat <- do.call(control[["get_prev.FUN"]], list(dat, at))
+ }
+
+
+ ## Track progress
+ verbose.icm(dat, type = "progress", s, at)
+ }
+
+ # Set output
+ if (s == 1) {
+ out <- saveout.seiqhrf.icm(dat, s)
+ } else {
+ out <- saveout.seiqhrf.icm(dat, s, out)
+ }
+
+ } # Simulation loop end
+
+ class(out) <- "icm"
+
+ } # end of single core execution
+
+ if (ncores > 1) {
+ doParallel::registerDoParallel(ncores)
+
+ sout <- foreach(s = 1:nsims) %dopar% {
+
+ control$nsims <- 1
+ control$currsim <- s
+
+ ## Initialization module
+ if (!is.null(control[["initialize.FUN"]])) {
+ dat <- do.call(control[["initialize.FUN"]], list(param, init, control))
+ }
+
+ # Timestep loop
+ for (at in 2:control$nsteps) {
+
+ ## User Modules
+ um <- control$user.mods
+ if (length(um) > 0) {
+ for (i in 1:length(um)) {
+ dat <- do.call(control[[um[i]]], list(dat, at))
+ }
+ }
+
+ ## Infection
+ if (!is.null(control[["infection.FUN"]])) {
+ dat <- do.call(control[["infection.FUN"]], list(dat, at))
+ }
+
+
+ ## Recovery
+ if (!is.null(control[["recovery.FUN"]])) {
+ dat <- do.call(control[["recovery.FUN"]], list(dat, at))
+ }
+
+
+ ## Departure Module
+ if (!is.null(control[["departures.FUN"]])) {
+ dat <- do.call(control[["departures.FUN"]], list(dat, at))
+ }
+
+
+ ## Arrival Module
+ if (!is.null(control[["arrivals.FUN"]])) {
+ dat <- do.call(control[["arrivals.FUN"]], list(dat, at))
+ }
+
+
+ ## Outputs
+ if (!is.null(control[["get_prev.FUN"]])) {
+ dat <- do.call(control[["get_prev.FUN"]], list(dat, at))
+ }
+
+
+ ## Track progress
+ verbose.icm(dat, type = "progress", s, at)
+ }
+
+ # Set output
+ out <- saveout.seiqhrf.icm(dat, s=1)
+ class(out) <- "icm"
+ return(out)
+
+ }
+
+ # aggregate results collected from each thread
+ collected_times <- list()
+
+ # collect the times from sout then delete them
+ for (i in 1:length(sout)) {
+ collected_times[[paste0("sim", i)]] <- sout[[i]]$times$sim1
+ sout[[i]]$times <- NULL
+ }
+
+ # merge $epi structures
+ merged.out <- sout[[1]]
+ for (i in 2:length(sout)) {
+ merged.out <- merge.seiqhrf.icm(merged.out, sout[[i]], param.error = FALSE)
+ }
+ out <- merged.out
+
+ # add the collected timing data
+ out$times <- collected_times
+
+ class(out) <- "icm"
+ } # end of parallel execution
+
+ return(out)
+}
+
+
+#' Initialize status of ICM
+#'
+#' Function to get the status of the initialized icm
+#'
+#' @param dat Object containing all data
+#'
+#' @return Updated dat
+init_status.icm <- function(dat) {
+
+ # Variables ---------------------------------------------------------------
+ type <- dat$control$type
+ group <- dat$attr$group
+ nGroups <- dat$param$groups
+
+ nG1 <- sum(group == 1)
+ nG2 <- sum(group == 2)
+
+ e.num <- dat$init$e.num
+ i.num <- dat$init$i.num
+ q.num <- dat$init$q.num
+ h.num <- dat$init$h.num
+ r.num <- dat$init$r.num
+ f.num <- dat$init$f.num
+ e.num.g2 <- dat$init$e.num.g2
+ i.num.g2 <- dat$init$i.num.g2
+ q.num.g2 <- dat$init$q.num.g2
+ h.num.g2 <- dat$init$h.num.g2
+ r.num.g2 <- dat$init$r.num.g2
+ f.num.g2 <- dat$init$f.num.g2
+
+ # Status ------------------------------------------------------------------
+ status <- rep("s", nG1 + nG2)
+ status[sample(which(group == 1), size = i.num)] <- "i"
+ if (nGroups == 2) {
+ status[sample(which(group == 2), size = i.num.g2)] <- "i"
+ }
+ if (type %in% c("SIR", "SEIR", "SEIQHR", "SEIQHRF")) {
+ status[sample(which(group == 1 & status == "s"), size = r.num)] <- "r"
+ if (nGroups == 2) {
+ status[sample(which(group == 2 & status == "s"), size = r.num.g2)] <- "r"
+ }
+ }
+ if (type %in% c("SEIR", "SEIQHR", "SEIQHRF")) {
+ status[sample(which(group == 1 & status == "s"), size = e.num)] <- "e"
+ if (nGroups == 2) {
+ status[sample(which(group == 2 & status == "s"), size = e.num.g2)] <- "e"
+ }
+ }
+ if (type %in% c("SEIQHR", "SEIQHRF")) {
+ status[sample(which(group == 1 & status == "s"), size = q.num)] <- "q"
+ if (nGroups == 2) {
+ status[sample(which(group == 2 & status == "s"), size = q.num.g2)] <- "q"
+ }
+ status[sample(which(group == 1 & status == "s"), size = h.num)] <- "h"
+ if (nGroups == 2) {
+ status[sample(which(group == 2 & status == "s"), size = h.num.g2)] <- "h"
+ }
+ }
+ if (type %in% c("SEIQHRF")) {
+ status[sample(which(group == 1 & status == "s"), size = f.num)] <- "f"
+ if (nGroups == 2) {
+ status[sample(which(group == 2 & status == "s"), size = f.num.g2)] <- "f"
+ }
+ }
+
+ dat$attr$status <- status
+
+
+ # Exposure Time ----------------------------------------------------------
+ idsExp <- which(status == "e")
+ expTime <- rep(NA, length(status))
+ # leave exposure time uninitialised for now, and
+ # just set to NA at start.
+ dat$attr$expTime <- expTime
+
+ # Infection Time ----------------------------------------------------------
+ idsInf <- which(status == "i")
+ infTime <- rep(NA, length(status))
+ dat$attr$infTime <- infTime # overwritten below
+
+ # Recovery Time ----------------------------------------------------------
+ idsRecov <- which(status == "r")
+ recovTime <- rep(NA, length(status))
+ dat$attr$recovTime <- recovTime
+
+ # Need for Hospitalisation Time ----------------------------------------------------------
+ idsHosp <- which(status == "h")
+ hospTime <- rep(NA, length(status))
+ dat$attr$hospTime <- hospTime
+
+ # Quarantine Time ----------------------------------------------------------
+ idsQuar <- which(status == "q")
+ quarTime <- rep(NA, length(status))
+ dat$attr$quarTime <- quarTime
+
+ # Hospital-need cessation Time ----------------------------------------------------------
+ dischTime <- rep(NA, length(status))
+ dat$attr$dischTime <- dischTime
+
+ # Case-fatality Time ----------------------------------------------------------
+ fatTime <- rep(NA, length(status))
+ dat$attr$fatTime <- fatTime
+
+ # If vital=TRUE, infTime is a uniform draw over the duration of infection
+ # note the initial infections may have negative infTime!
+ if (FALSE) {
+ # not sure what the following section is trying to do, but it
+ # mucks up the gamma-distributed incumabtion periods, so set
+ # infTime for initial infected people to t=1 instead
+ if (dat$param$vital == TRUE && dat$param$di.rate > 0) {
+ infTime[idsInf] <- -rgeom(n = length(idsInf), prob = dat$param$di.rate) + 2
+ } else {
+ if (dat$control$type == "SI" || dat$param$rec.rate == 0) {
+ # infTime a uniform draw over the number of sim time steps
+ infTime[idsInf] <- ssample(1:(-dat$control$nsteps + 2),
+ length(idsInf), replace = TRUE)
+ } else {
+ if (nGroups == 1) {
+ infTime[idsInf] <- ssample(1:(-round(1 / dat$param$rec.rate) + 2),
+ length(idsInf), replace = TRUE)
+ }
+ if (nGroups == 2) {
+ infG1 <- which(status == "i" & group == 1)
+ infTime[infG1] <- ssample(1:(-round(1 / dat$param$rec.rate) + 2),
+ length(infG1), replace = TRUE)
+ infG2 <- which(status == "i" & group == 2)
+ infTime[infG2] <- ssample(1:(-round(1 / dat$param$rec.rate.g2) + 2),
+ length(infG2), replace = TRUE)
+ }
+ }
+ }
+ }
+ infTime[idsInf] <- 1
+ dat$attr$infTime <- infTime
+
+ return(dat)
+}
+
+
+#' Check control ICM
+#'
+#' Function to check the control icm
+#'
+#' @param type Model type
+#'
+#' @return Control icm
+control.icm <- function(type, nsteps, nsims = 1,
+ rec.rand = TRUE, quar.rand = TRUE, hosp.rand = TRUE, disch.rand = TRUE,
+ fat.rand = TRUE, a.rand = TRUE, d.rand = TRUE, initialize.FUN = initialize.FUN,
+ infection.FUN = infection.FUN, recovery.FUN = recovery.FUN,
+ departures.FUN = departures.FUN, arrivals.FUN = arrivals.FUN,
+ get_prev.FUN = get_prev.FUN, verbose = FALSE,
+ verbose.int = 0, skip.check = FALSE, ncores=1, ...) {
+
+ # Get arguments
+ p <- list()
+ formal.args <- formals(sys.function())
+ formal.args[["..."]] <- NULL
+ for (arg in names(formal.args)) {
+ if (as.logical(mget(arg) != "")) {
+ p[arg] <- list(get(arg))
+ }
+ }
+ dot.args <- list(...)
+ names.dot.args <- names(dot.args)
+ if (length(dot.args) > 0) {
+ for (i in 1:length(dot.args)) {
+ p[[names.dot.args[i]]] <- dot.args[[i]]
+ }
+ }
+
+ if ("births.FUN" %in% names(dot.args)) {
+ p$arrivals.FUN <- dot.args$births.FUN
+ p$births.FUN <- dot.args$births.FUN <- NULL
+ message("EpiModel 1.7.0 onward renamed the birth function births.FUN to arrivals.FUN. See documentation for details.")
+ }
+ if ("deaths.FUN" %in% names(dot.args)) {
+ p$departures.FUN <- dot.args$deaths.FUN
+ p$deaths.FUN <- dot.args$deaths.FUN <- NULL
+ message("EpiModel 1.7.0 onward renamed the death function deaths.FUN to departures.FUN. See documentation for details.")
+ }
+
+
+ ## Module classification
+ p$bi.mods <- grep(".FUN", names(formal.args), value = TRUE)
+ p$user.mods <- grep(".FUN", names(dot.args), value = TRUE)
+
+
+ ## Defaults and checks
+ if (is.null(p$type) | !(p$type %in% c("SI", "SIS", "SIR", "SEIR", "SEIQHR", "SEIQHRF"))) {
+ stop("Specify type as \"SI\", \"SIS\", \"SIR\", \"SEIR\", \"SEIQHR\", or \"SEIQHRF\" ", call. = FALSE)
+ }
+ if (is.null(p$nsteps)) {
+ stop("Specify nsteps", call. = FALSE)
+ }
+
+
+ ## Output
+ class(p) <- c("control.icm", "list")
+ return(p)
+}
+
+
+#' Merge icm
+#'
+#' Function to merge icms
+#'
+#' @param x First to merge
+#' @param y Second to merge
+#'
+#' @return Combined
+merge.seiqhrf.icm <- function(x, y, ...) {
+
+ ## Check structure
+ if (length(x) != length(y) || !identical(names(x), names(y))) {
+ stop("x and y have different structure")
+ }
+ if (x$control$nsims > 1 & y$control$nsims > 1 &
+ !identical(sapply(x, class), sapply(y, class))) {
+ stop("x and y have different structure")
+ }
+
+ ## Check params
+ check1 <- identical(x$param, y$param)
+ check2 <- identical(x$control[-which(names(x$control) %in% c("nsims", "currsim"))],
+ y$control[-which(names(y$control) %in% c("nsims", "currsim"))])
+
+ if (check1 == FALSE) {
+ stop("x and y have different parameters")
+ }
+ if (check2 == FALSE) {
+ stop("x and y have different controls")
+ }
+
+
+ z <- x
+ new.range <- (x$control$nsims + 1):(x$control$nsims + y$control$nsims)
+
+ # Merge data
+ for (i in 1:length(x$epi)) {
+ if (x$control$nsims == 1) {
+ x$epi[[i]] <- data.frame(x$epi[[i]])
+ }
+ if (y$control$nsims == 1) {
+ y$epi[[i]] <- data.frame(y$epi[[i]])
+ }
+ z$epi[[i]] <- cbind(x$epi[[i]], y$epi[[i]])
+ names(z$epi[[i]])[new.range] <- paste0("sim", new.range)
+ }
+
+ z$control$nsims <- max(new.range)
+
+ return(z)
+}
+
+
+
+#' Progress icm
+#'
+#' Function to get progress of icms
+#'
+#' @param dat Object containing all data
+#' @param at ?
+#'
+#' @return progress
+progress.seiqhrf.icm <- function(dat, at) {
+
+ #print(at)
+ #print(dat$control$type)
+ #print("-------")
+
+ # Conditions --------------------------------------------------------------
+ if (!(dat$control$type %in% c("SIR", "SIS", "SEIR", "SEIQHR", "SEIQHRF"))) {
+ return(dat)
+ }
+
+
+ # Variables ---------------------------------------------------------------
+ active <- dat$attr$active
+ status <- dat$attr$status
+
+ groups <- dat$param$groups
+ group <- dat$attr$group
+
+ type <- dat$control$type
+ recovState <- ifelse(type %in% c("SIR", "SEIR", "SEIQHR", "SEIQHRF"), "r", "s")
+ progState <- "i"
+ quarState <- "q"
+ hospState <- "h"
+ fatState <- "f"
+
+ # --- progress from exposed to infectious ----
+ prog.rand <- dat$control$prog.rand
+ prog.rate <- dat$param$prog.rate
+ prog.rate.g2 <- dat$param$prog.rate.g2
+ prog.dist.scale <- dat$param$prog.dist.scale
+ prog.dist.shape <- dat$param$prog.dist.shape
+ prog.dist.scale.g2 <- dat$param$prog.dist.scale.g2
+ prog.dist.shape.g2 <- dat$param$prog.dist.shape.g2
+
+ nProg <- nProgG2 <- 0
+ idsElig <- which(active == 1 & status == "e")
+ nElig <- length(idsElig)
+
+ if (nElig > 0) {
+
+ gElig <- group[idsElig]
+ rates <- c(prog.rate, prog.rate.g2)
+ ratesElig <- rates[gElig]
+
+ if (prog.rand == TRUE) {
+ vecProg <- which(rbinom(nElig, 1, ratesElig) == 1)
+ if (length(vecProg) > 0) {
+ idsProg <- idsElig[vecProg]
+ nProg <- sum(group[idsProg] == 1)
+ nProgG2 <- sum(group[idsProg] == 2)
+ status[idsProg] <- progState
+ dat$attr$infTime[idsProg] <- at
+ }
+ } else {
+ vecTimeSinceExp <- at - dat$attr$expTime[idsElig]
+ gammaRatesElig <- pweibull(vecTimeSinceExp, prog.dist.shape, scale=prog.dist.scale)
+ nProg <- round(sum(gammaRatesElig[gElig == 1], na.rm=TRUE))
+ if (nProg > 0) {
+ ids2bProg <- ssample(idsElig[gElig == 1],
+ nProg, prob = gammaRatesElig[gElig == 1])
+ status[ids2bProg] <- progState
+ dat$attr$infTime[ids2bProg] <- at
+ # debug
+ if (FALSE & at <= 30) {
+ print(paste("at:", at))
+ print("idsElig:")
+ print(idsElig[gElig == 1])
+ print("vecTimeSinceExp:")
+ print(vecTimeSinceExp[gElig == 1])
+ print("gammaRatesElig:")
+ print(gammaRatesElig)
+ print(paste("nProg:",nProg))
+ print(paste("sum of elig rates:", round(sum(gammaRatesElig[gElig == 1]))))
+ print(paste("sum(gElig == 1):", sum(gElig == 1)))
+ print("ids progressed:")
+ print(ids2bProg)
+ print("probs of ids to be progressed:")
+ print(gammaRatesElig[which(idsElig %in% ids2bProg)])
+ print("days since exposed of ids to be progressed:")
+ print(vecTimeSinceExp[which(idsElig %in% ids2bProg)])
+ print("------")
+ }
+ }
+ if (groups == 2) {
+ nProgG2 <- round(sum(gammaRatesElig[gElig == 2], na.rm=TRUE))
+ if (nProgG2 > 0) {
+ ids2bProgG2 <- ssample(idsElig[gElig == 2],
+ nProgG2, prob = gammaRatesElig[gElig == 2])
+ status[ids2bProgG2] <- progState
+ dat$attr$infTime[ids2bProgG2] <- at
+ }
+ }
+ }
+ }
+ dat$attr$status <- status
+
+ if (type %in% c("SEIQHR", "SEIQHRF")) {
+ # ----- quarantine -------
+ quar.rand <- dat$control$quar.rand
+ quar.rate <- dat$param$quar.rate
+ quar.rate.g2 <- dat$param$quar.rate.g2
+
+ nQuar <- nQuarG2 <- 0
+ idsElig <- which(active == 1 & status == "i")
+ nElig <- length(idsElig)
+
+ if (nElig > 0) {
+
+ gElig <- group[idsElig]
+ rates <- c(quar.rate, quar.rate.g2)
+
+ if (length(quar.rate) > 1) {
+ qrate <- quar.rate[at]
+ } else {
+ qrate <- quar.rate
+ }
+ if (length(quar.rate.g2) > 1) {
+ qrate.g2 <- quar.rate.g2[at]
+ } else {
+ qrate.g2 <- quar.rate.g2
+ }
+ rates <- c(qrate, qrate.g2)
+ ratesElig <- rates[gElig]
+ if (quar.rand == TRUE) {
+ vecQuar <- which(rbinom(nElig, 1, ratesElig) == 1)
+ if (length(vecQuar) > 0) {
+ idsQuar <- idsElig[vecQuar]
+ nQuar <- sum(group[idsQuar] == 1)
+ nQuarG2 <- sum(group[idsQuar] == 2)
+ status[idsQuar] <- quarState
+ dat$attr$quarTime[idsQuar] <- at
+ }
+ } else {
+ nQuar <- min(round(sum(ratesElig[gElig == 1])), sum(gElig == 1))
+ idsQuar <- ssample(idsElig[gElig == 1], nQuar)
+ status[idsQuar] <- quarState
+ dat$attr$quarTime[idsQuar] <- at
+ if (groups == 2) {
+ nQuarG2 <- min(round(sum(ratesElig[gElig == 2])), sum(gElig == 2))
+ idsQuarG2 <- ssample(idsElig[gElig == 2], nQuarG2)
+ status[idsQuarG2] <- quarState
+ dat$attr$quarTime[idsQuarG2] <- at
+ }
+ }
+ }
+ dat$attr$status <- status
+
+ # ----- need to be hospitalised -------
+ hosp.rand <- dat$control$hosp.rand
+ hosp.rate <- dat$param$hosp.rate
+ hosp.rate.g2 <- dat$param$hosp.rate.g2
+
+ nHosp <- nHospG2 <- 0
+ idsElig <- which(active == 1 & (status == "i" | status == "q"))
+ nElig <- length(idsElig)
+ idsHosp <- numeric(0)
+
+ if (nElig > 0) {
+
+ gElig <- group[idsElig]
+ rates <- c(hosp.rate, hosp.rate.g2)
+ ratesElig <- rates[gElig]
+
+ if (hosp.rand == TRUE) {
+ vecHosp <- which(rbinom(nElig, 1, ratesElig) == 1)
+ if (length(vecHosp) > 0) {
+ idsHosp <- idsElig[vecHosp]
+ nHosp <- sum(group[idsHosp] == 1)
+ nHospG2 <- sum(group[idsHosp] == 2)
+ status[idsHosp] <- hospState
+ }
+ } else {
+ nHosp <- min(round(sum(ratesElig[gElig == 1])), sum(gElig == 1))
+ idsHosp <- ssample(idsElig[gElig == 1], nHosp)
+ status[idsHosp] <- hospState
+ if (groups == 2) {
+ nHospG2 <- min(round(sum(ratesElig[gElig == 2])), sum(gElig == 2))
+ idsHospG2 <- ssample(idsElig[gElig == 2], nHospG2)
+ status[idsHospG2] <- hospState
+ idsHosp <- c(idsHosp, idsHospG2)
+ }
+ }
+ }
+ dat$attr$status <- status
+ dat$attr$hospTime[idsHosp] <- at
+
+ # ----- discharge from need to be hospitalised -------
+ disch.rand <- dat$control$disch.rand
+ disch.rate <- dat$param$disch.rate
+ disch.rate.g2 <- dat$param$disch.rate.g2
+
+ nDisch <- nDischG2 <- 0
+ idsElig <- which(active == 1 & status == "h")
+ nElig <- length(idsElig)
+ idsDisch <- numeric(0)
+
+ if (nElig > 0) {
+
+ gElig <- group[idsElig]
+ rates <- c(disch.rate, disch.rate.g2)
+
+ if (length(disch.rate) > 1) {
+ dcrate <- disch.rate[at]
+ } else {
+ dcrate <- disch.rate
+ }
+ if (length(disch.rate.g2) > 1) {
+ dcrate.g2 <- disch.rate.g2[at]
+ } else {
+ dcrate.g2 <- disch.rate.g2
+ }
+
+ rates <- c(dcrate, dcrate.g2)
+ ratesElig <- rates[gElig]
+
+ if (disch.rand == TRUE) {
+ vecDisch <- which(rbinom(nElig, 1, ratesElig) == 1)
+ if (length(vecDisch) > 0) {
+ idsDisch <- idsElig[vecDisch]
+ nDisch <- sum(group[idsDisch] == 1)
+ nDischG2 <- sum(group[idsDisch] == 2)
+ status[idsDisch] <- recovState
+ }
+ } else {
+ nDisch <- min(round(sum(ratesElig[gElig == 1])), sum(gElig == 1))
+ idsDisch <- ssample(idsElig[gElig == 1], nDisch)
+ status[idsDisch] <- recovState
+ if (groups == 2) {
+ nDischG2 <- min(round(sum(ratesElig[gElig == 2])), sum(gElig == 2))
+ idsDischG2 <- ssample(idsElig[gElig == 2], nDischG2)
+ status[idsDischG2] <- recovState
+ idsDisch <- c(idsDisch, idsDischG2)
+ }
+ }
+ }
+ dat$attr$status <- status
+ dat$attr$dischTime[idsDisch] <- at
+ }
+
+ # ----- recover -------
+ rec.rand <- dat$control$rec.rand
+ rec.rate <- dat$param$rec.rate
+ rec.rate.g2 <- dat$param$rec.rate.g2
+ rec.dist.scale <- dat$param$rec.dist.scale
+ rec.dist.shape <- dat$param$rec.dist.shape
+ rec.dist.scale.g2 <- dat$param$rec.dist.scale.g2
+ rec.dist.shape.g2 <- dat$param$rec.dist.shape.g2
+
+ nRecov <- nRecovG2 <- 0
+ idsElig <- which(active == 1 & (status == "i" | status == "q" | status == "h"))
+ nElig <- length(idsElig)
+ idsRecov <- numeric(0)
+
+ if (nElig > 0) {
+
+ gElig <- group[idsElig]
+ rates <- c(rec.rate, rec.rate.g2)
+ ratesElig <- rates[gElig]
+
+ if (rec.rand == TRUE) {
+ vecRecov <- which(rbinom(nElig, 1, ratesElig) == 1)
+ if (length(vecRecov) > 0) {
+ idsRecov <- idsElig[vecRecov]
+ nRecov <- sum(group[idsRecov] == 1)
+ nRecovG2 <- sum(group[idsRecov] == 2)
+ status[idsRecov] <- recovState
+ }
+ } else {
+ vecTimeSinceExp <- at - dat$attr$expTime[idsElig]
+ vecTimeSinceExp[is.na(vecTimeSinceExp)] <- 0
+ gammaRatesElig <- pweibull(vecTimeSinceExp, rec.dist.shape, scale=rec.dist.scale)
+ nRecov <- round(sum(gammaRatesElig[gElig == 1], na.rm=TRUE))
+ if (nRecov > 0) {
+ idsRecov <- ssample(idsElig[gElig == 1],
+ nRecov, prob = gammaRatesElig[gElig == 1])
+ status[idsRecov] <- recovState
+ # debug
+ if (FALSE & at <= 30) {
+ print(paste("at:", at))
+ print("idsElig:")
+ print(idsElig[gElig == 1])
+ print("vecTimeSinceExp:")
+ print(vecTimeSinceExp[gElig == 1])
+ print("gammaRatesElig:")
+ print(gammaRatesElig)
+ print(paste("nRecov:",nRecov))
+ print(paste("sum of elig rates:", round(sum(gammaRatesElig[gElig == 1]))))
+ print(paste("sum(gElig == 1):", sum(gElig == 1)))
+ print("ids recovered:")
+ print(idsRecov)
+ print("probs of ids to be progressed:")
+ print(gammaRatesElig[which(idsElig %in% idsRecov)])
+ print("days since exposed of ids to be Recovered:")
+ print(vecTimeSinceExp[which(idsElig %in% idsRecov)])
+ print("------")
+ }
+
+ }
+ if (groups == 2) {
+ nRecovG2 <- round(sum(gammaRatesElig[gElig == 2], na.rm=TRUE))
+ if (nRecovG2 > 0) {
+ idsRecovG2 <- ssample(idsElig[gElig == 2],
+ nRecovG2, prob = gammaRatesElig[gElig == 2])
+ status[idsRecovG2] <- recovState
+ idsRecov <- c(idsRecov, idsRecovG2)
+ }
+ }
+ }
+ }
+ dat$attr$status <- status
+ dat$attr$recovTime[idsRecov] <- at
+
+ fatEnable <- TRUE
+ if (fatEnable & type %in% c("SEIQHRF")) {
+ # ----- case fatality -------
+ fat.rand <- dat$control$fat.rand
+ fat.rate.base <- dat$param$fat.rate.base
+ fat.rate.base.g2 <- dat$param$fat.rate.base.g2
+ fat.rate.base.g2 <- ifelse(is.null(fat.rate.base.g2),
+ 0, fat.rate.base.g2)
+ fat.rate.overcap <- dat$param$fat.rate.overcap
+ fat.rate.overcap.g2 <- dat$param$fat.rate.overcap.g2
+ fat.rate.overcap.g2 <- ifelse(is.null(fat.rate.overcap.g2),
+ 0, fat.rate.overcap.g2)
+ hosp.cap <- dat$param$hosp.cap
+ fat.tcoeff <- dat$param$fat.tcoeff
+
+ nFat <- nFatG2 <- 0
+ idsElig <- which(active == 1 & status =="h")
+ nElig <- length(idsElig)
+
+ if (nElig > 0) {
+ gElig <- group[idsElig]
+ timeInHospElig <- at - dat$attr$hospTime[idsElig]
+ rates <- c(fat.rate.base, fat.rate.base.g2)
+ h.num.yesterday <- 0
+ if (!is.null(dat$epi$h.num[at - 1])) {
+ h.num.yesterday <- dat$epi$h.num[at - 1]
+ if (h.num.yesterday > hosp.cap) {
+ blended.rate <- ((hosp.cap * fat.rate.base) +
+ ((h.num.yesterday - hosp.cap) * fat.rate.overcap)) /
+ h.num.yesterday
+ blended.rate.g2 <- ((hosp.cap * fat.rate.base.g2) +
+ ((h.num.yesterday - hosp.cap) * fat.rate.overcap.g2)) /
+ h.num.yesterday
+ rates <- c(blended.rate, blended.rate.g2)
+ }
+ }
+ ratesElig <- rates[gElig]
+ ratesElig <- ratesElig + timeInHospElig*fat.tcoeff*ratesElig
+
+ if (fat.rand == TRUE) {
+ vecFat <- which(rbinom(nElig, 1, ratesElig) == 1)
+ if (length(vecFat) > 0) {
+ idsFat <- idsElig[vecFat]
+ nFat <- sum(group[idsFat] == 1)
+ nFatG2 <- sum(group[idsFat] == 2)
+ status[idsFat] <- fatState
+ dat$attr$fatTime[idsFat] <- at
+ }
+ } else {
+ nFat <- min(round(sum(ratesElig[gElig == 1])), sum(gElig == 1))
+ idsFat <- ssample(idsElig[gElig == 1], nFat)
+ status[idsFat] <- fatState
+ dat$attr$fatTime[idsFat] <- at
+ if (groups == 2) {
+ nFatG2 <- min(round(sum(ratesElig[gElig == 2])), sum(gElig == 2))
+ idsFatG2 <- ssample(idsElig[gElig == 2], nFatG2)
+ status[idsFatG2] <- fatState
+ dat$attr$fatTime[idsFatG2] <- at
+ }
+ }
+ }
+ dat$attr$status <- status
+ }
+
+ # Output ------------------------------------------------------------------
+ outName_a <- ifelse(type %in% c("SIR", "SEIR"), "ir.flow", "is.flow")
+ outName_a[2] <- paste0(outName_a, ".g2")
+ if (type %in% c("SEIR", "SEIQHR", "SEIQHRF")) {
+ outName_b <- "ei.flow"
+ outName_b[2] <- paste0(outName_b, ".g2")
+ }
+ if (type %in% c("SEIQHR", "SEIQHRF")) {
+ outName_c <- "iq.flow"
+ outName_c[2] <- paste0(outName_c, ".g2")
+ outName_d <- "iq2h.flow"
+ outName_d[2] <- paste0(outName_d, ".g2")
+ }
+ if (type %in% c("SEIQHRF")) {
+ outName_e <- "hf.flow"
+ outName_e[2] <- paste0(outName_e, ".g2")
+ }
+ ## Summary statistics
+ if (at == 2) {
+ dat$epi[[outName_a[1]]] <- c(0, nRecov)
+ if (type %in% c("SEIR", "SEIQHR")) {
+ dat$epi[[outName_b[1]]] <- c(0, nProg)
+ }
+ if (type %in% c("SEIQHR", "SEIQHRF")) {
+ dat$epi[[outName_c[1]]] <- c(0, nQuar)
+ dat$epi[[outName_d[1]]] <- c(0, nHosp)
+ }
+ if (fatEnable & type %in% c("SEIQHRF")) {
+ dat$epi[[outName_e[1]]] <- c(0, nFat)
+ }
+ } else {
+ dat$epi[[outName_a[1]]][at] <- nRecov
+ if (type %in% c("SEIR", "SEIQHR")) {
+ dat$epi[[outName_b[1]]][at] <- nProg
+ }
+ if (type %in% c("SEIQHR", "SEIQHRF")) {
+ dat$epi[[outName_c[1]]][at] <- nQuar
+ dat$epi[[outName_d[1]]][at] <- nHosp
+ }
+ if (fatEnable & type %in% c("SEIQHRF")) {
+ dat$epi[[outName_e[1]]][at] <- nFat
+ }
+ }
+ if (groups == 2) {
+ if (at == 2) {
+ dat$epi[[outName_a[2]]] <- c(0, nRecovG2)
+ if (type %in% c("SEIR", "SEIQHR", "SEIQHRF")) {
+ dat$epi[[outName_b[2]]] <- c(0, nProgG2)
+ }
+ if (type %in% c("SEIQHR", "SEIQHRF")) {
+ dat$epi[[outName_c[2]]] <- c(0, nQuarG2)
+ dat$epi[[outName_d[2]]] <- c(0, nHospG2)
+ }
+ if (type %in% c("SEIQHRF")) {
+ dat$epi[[outName_e[2]]] <- c(0, nFatG2)
+ }
+ } else {
+ dat$epi[[outName_a[2]]][at] <- nRecovG2
+ if (type %in% c("SEIR", "SEIQHR", "SEIQHRF")) {
+ dat$epi[[outName_b[2]]][at] <- nProgG2
+ }
+ if (type %in% c("SEIQHR", "SEIQHRF")) {
+ dat$epi[[outName_c[2]]][at] <- nQuarG2
+ dat$epi[[outName_d[2]]][at] <- nHospG2
+ }
+ if (type %in% c("SEIQHRF")) {
+ dat$epi[[outName_e[2]]][at] <- nFatG2
+ }
+ }
+ }
+
+ return(dat)
+}
+
+#' Save icm
+#'
+#' Function to save icm
+#'
+#' @param dat Object containing all data
+#' @param at ?
+#' @param out ?
+#'
+#' @return out to save
+#'
+saveout.seiqhrf.icm <- function(dat, s, out = NULL) {
+
+ alist2df <- function(dat,s) {
+ alist <- list()
+ alist$expTime <- dat$attr$expTime
+ alist$infTime <- dat$attr$infTime
+ alist$quarTime <- dat$attr$quarTime
+ alist$recovTime <- dat$attr$recovTime
+ alist$hospTime <- dat$attr$hospTime
+ alist$dischTime <- dat$attr$dischTime
+ alist$fatTime <- dat$attr$fatTime
+ alist <- lapply(alist, `length<-`, max(lengths(alist)))
+ return(data.frame(alist))
+ }
+
+ if (s == 1) {
+ out <- list()
+ out$param <- dat$param
+ out$control <- dat$control
+ out$epi <- list()
+ for (j in 1:length(dat$epi)) {
+ out$epi[[names(dat$epi)[j]]] <- data.frame(dat$epi[j])
+ }
+ } else {
+ for (j in 1:length(dat$epi)) {
+ out$epi[[names(dat$epi)[j]]][, s] <- data.frame(dat$epi[j])
+ }
+ }
+
+ # capture transition times from attribs
+ if (dat$control$type %in% c("SEIQHR", "SEIQHRF")) {
+ out$times[[paste("sim",s,sep="")]] <- alist2df(dat,s)
+ }
+
+ ## Processing for final run
+ if (s == dat$getParam(control, 'nsims')) {
+
+ # Remove functions from control list
+ ftodel <- grep(".FUN", names(out$control), value = TRUE)
+ out$control[ftodel] <- NULL
+
+ # Set column names for varying list elements
+ for (i in as.vector(which(lapply(out$epi, class) == "data.frame"))) {
+ colnames(out$epi[[i]]) <- paste0("sim", 1:dat$getParam(control, 'nsims'))
+ }
+
+ }
+
+ return(out)
+}
+
+
\ No newline at end of file
diff --git a/R/sirPlus-package.R b/R/sirPlus-package.R
new file mode 100644
index 0000000..49512df
--- /dev/null
+++ b/R/sirPlus-package.R
@@ -0,0 +1,9 @@
+#' sirPlus
+#'
+#' \pkg{sirPlus} is a package for modeling COVID-19 spread, focusing on rates
+#' of hospitalization, which should be useful for hospitals.
+#'
+#'
+#' @name sirPlus
+#' @docType package
+NULL
diff --git a/R/utils.R b/R/utils.R
index 155ea80..30886a1 100644
--- a/R/utils.R
+++ b/R/utils.R
@@ -1,38 +1,3 @@
-#' @importFrom utils globalVariables
-utils::globalVariables(c("pval_nominal", 'perc', 'ex_pairs', 'ex_means'))
-
-#' Logistic function
-#'
-#' Implementation of the logistic function
-#'
-#' @param x value to apply the function to.
-#' @param x0 midpoint parameter. Gives the centre of the function.
-#' @param k shape parameter. Gives the slope of the function.
-#'
-#' @return Value of logistic function with given parameters
-logistic <- function(x, x0, k) {
- 1 / (1 + exp(-k * (x - x0)))
-}
-
-#' Bind rows (matched)
-#'
-#' Bind the rows of two data frames, keeping only the columns that are common
-#' to both.
-#'
-#' @param df1 first data.frame to bind.
-#' @param df2 second data.frame to bind.
-#'
-#' @return data.frame containing rows from \code{df1} and \code{df2} but only
-#' common columns.
-rbindMatched <- function(df1, df2) {
- common.names <- intersect(names(df1), names(df2))
- if (length(common.names) < 2) {
- stop("There must be at least two columns in common")
- }
- combined <- rbind(df1[, common.names], df2[, common.names])
-
- return(combined)
-}
#' Bring items forward
#'
@@ -61,39 +26,28 @@ bringItemsForward <- function(ll, items) {
return(ll)
}
-#' Winsorize vector
-#'
-#' Set outliers in a numeric vector to a specified percentile.
-#'
-#' @param x Numeric vector to winsorize
-#' @param q Percentile to set from each end
-#'
-#' @return Winsorized numeric vector
-winsorize <- function(x, q) {
-
- checkmate::check_numeric(x, any.missing = FALSE)
- checkmate::check_number(q, lower = 0, upper = 1)
-
- lohi <- stats::quantile(x, c(q, 1 - q), na.rm = TRUE)
-
- if (diff(lohi) < 0) { lohi <- rev(lohi) }
-
- x[!is.na(x) & x < lohi[1]] <- lohi[1]
- x[!is.na(x) & x > lohi[2]] <- lohi[2]
-
- return(x)
-}
-
-#' Calculate coefficient of variation
-#'
-#' Implementation of the coefficient of variation
+#' cum_discr_si
#'
-#' @param x vector of values.
+#' cum_discr_si
#'
+#' @param vecTimeSinceExp ?
+#' @param scale ?
+#' @param shape ?
+#'
#' @return Value of coefficient of variation for vector
#' @importFrom stats sd
-co.var <- function(x) {
- sd(x)/mean(x)
-}
+cum_discr_si <- function(vecTimeSinceExp, scale, shape) {
+ vlen <- length(vecTimeSinceExp)
+ if (vlen > 0) {
+ probVec <- numeric(vlen)
+ for (p in 1:vlen) {
+ probVec[p] <- pweibull(vecTimeSinceExp[p], shape=shape, scale=scale)
+ }
+ } else {
+ probVec <- 0
+ }
+ return(probVec)
+}
+
diff --git a/man/arrivals.FUN.Rd b/man/arrivals.FUN.Rd
new file mode 100644
index 0000000..285952a
--- /dev/null
+++ b/man/arrivals.FUN.Rd
@@ -0,0 +1,23 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/epi-functions.R
+\name{arrivals.FUN}
+\alias{arrivals.FUN}
+\title{Arrivals function}
+\usage{
+arrivals.FUN(dat, at)
+}
+\arguments{
+\item{dat}{Input data needed.}
+
+\item{at}{Not sure???}
+}
+\value{
+Updated dat
+}
+\description{
+Function to handel background demographics for the SEIQHRF model.
+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..
+}
diff --git a/man/control.icm.Rd b/man/control.icm.Rd
new file mode 100644
index 0000000..92fe5aa
--- /dev/null
+++ b/man/control.icm.Rd
@@ -0,0 +1,39 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/sirPlus-functions.R
+\name{control.icm}
+\alias{control.icm}
+\title{Check control ICM}
+\usage{
+control.icm(
+ type,
+ nsteps,
+ nsims = 1,
+ rec.rand = TRUE,
+ quar.rand = TRUE,
+ hosp.rand = TRUE,
+ disch.rand = TRUE,
+ fat.rand = TRUE,
+ a.rand = TRUE,
+ d.rand = TRUE,
+ initialize.FUN = initialize.icm,
+ infection.FUN = infection.icm,
+ recovery.FUN = recovery.icm,
+ departures.FUN = departures.icm,
+ arrivals.FUN = arrivals.icm,
+ get_prev.FUN = get_prev.icm,
+ verbose = FALSE,
+ verbose.int = 0,
+ skip.check = FALSE,
+ ncores = 1,
+ ...
+)
+}
+\arguments{
+\item{type}{Model type}
+}
+\value{
+Control icm
+}
+\description{
+Function to check the control icm
+}
diff --git a/man/cum_discr_si.Rd b/man/cum_discr_si.Rd
new file mode 100644
index 0000000..18214e7
--- /dev/null
+++ b/man/cum_discr_si.Rd
@@ -0,0 +1,21 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/utils.R
+\name{cum_discr_si}
+\alias{cum_discr_si}
+\title{cum_discr_si}
+\usage{
+cum_discr_si(vecTimeSinceExp, scale, shape)
+}
+\arguments{
+\item{vecTimeSinceExp}{?}
+
+\item{scale}{?}
+
+\item{shape}{?}
+}
+\value{
+Value of coefficient of variation for vector
+}
+\description{
+cum_discr_si
+}
diff --git a/man/departures.FUN.Rd b/man/departures.FUN.Rd
new file mode 100644
index 0000000..b8aba6e
--- /dev/null
+++ b/man/departures.FUN.Rd
@@ -0,0 +1,20 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/epi-functions.R
+\name{departures.FUN}
+\alias{departures.FUN}
+\title{Departures function}
+\usage{
+departures.FUN(dat, at)
+}
+\arguments{
+\item{dat}{Input data needed.}
+
+\item{at}{Not sure???}
+}
+\value{
+Updated dat
+}
+\description{
+Function to handel background demographics for the SEIQHRF model.
+Specifically departures (deaths not due to the virus, and emigration).
+}
diff --git a/man/get_prev.FUN.Rd b/man/get_prev.FUN.Rd
new file mode 100644
index 0000000..eb0ad9f
--- /dev/null
+++ b/man/get_prev.FUN.Rd
@@ -0,0 +1,20 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/epi-functions.R
+\name{get_prev.FUN}
+\alias{get_prev.FUN}
+\title{Get previous function}
+\usage{
+get_prev.FUN(dat, at)
+}
+\arguments{
+\item{dat}{Input data needed.}
+
+\item{at}{Not sure???}
+}
+\value{
+Updated dat
+}
+\description{
+Utility function that collects prevalence and transition time data from each
+run and stores it away in the simulation result object.
+}
diff --git a/man/get_times.Rd b/man/get_times.Rd
new file mode 100644
index 0000000..2a43e33
--- /dev/null
+++ b/man/get_times.Rd
@@ -0,0 +1,17 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/results-parse.R
+\name{get_times}
+\alias{get_times}
+\title{Get times}
+\usage{
+get_times(simulate_results)
+}
+\arguments{
+\item{simulate_results}{results from `simulte()` function.}
+}
+\value{
+Dataframe of
+}
+\description{
+Function to extract timings and assemble results in a dataframe
+}
diff --git a/man/icm.seiqhrf.Rd b/man/icm.seiqhrf.Rd
new file mode 100644
index 0000000..4e1c29b
--- /dev/null
+++ b/man/icm.seiqhrf.Rd
@@ -0,0 +1,21 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/sirPlus-functions.R
+\name{icm.seiqhrf}
+\alias{icm.seiqhrf}
+\title{Infection function}
+\usage{
+icm.seiqhrf(param, init, control)
+}
+\arguments{
+\item{param}{ICM parameters.}
+
+\item{init}{Initial value parameters.}
+
+\item{control}{Control parameters.}
+}
+\value{
+Updated dat
+}
+\description{
+Function to implement infection processes for SEIQHRF model
+}
diff --git a/man/infection.FUN.Rd b/man/infection.FUN.Rd
new file mode 100644
index 0000000..2890047
--- /dev/null
+++ b/man/infection.FUN.Rd
@@ -0,0 +1,19 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/epi-functions.R
+\name{infection.FUN}
+\alias{infection.FUN}
+\title{Infection function}
+\usage{
+infection.FUN(dat, at)
+}
+\arguments{
+\item{dat}{Input data needed.}
+
+\item{at}{Not sure???}
+}
+\value{
+Updated dat
+}
+\description{
+Function to implement infection processes for SEIQHRF model
+}
diff --git a/man/infection.seiqhrf.icm.Rd b/man/infection.seiqhrf.icm.Rd
new file mode 100644
index 0000000..658db31
--- /dev/null
+++ b/man/infection.seiqhrf.icm.Rd
@@ -0,0 +1,19 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/sirPlus-functions.R
+\name{infection.seiqhrf.icm}
+\alias{infection.seiqhrf.icm}
+\title{Model infections in ICM}
+\usage{
+infection.seiqhrf.icm(dat, at)
+}
+\arguments{
+\item{dat}{Object containing all data}
+
+\item{at}{?}
+}
+\value{
+Updated dat
+}
+\description{
+Function to model infections in icm
+}
diff --git a/man/init_status.icm.Rd b/man/init_status.icm.Rd
new file mode 100644
index 0000000..1328bb6
--- /dev/null
+++ b/man/init_status.icm.Rd
@@ -0,0 +1,17 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/sirPlus-functions.R
+\name{init_status.icm}
+\alias{init_status.icm}
+\title{Initialize status of ICM}
+\usage{
+init_status.icm(dat)
+}
+\arguments{
+\item{dat}{Object containing all data}
+}
+\value{
+Updated dat
+}
+\description{
+Function to get the status of the initialized icm
+}
diff --git a/man/initialize.FUN.Rd b/man/initialize.FUN.Rd
new file mode 100644
index 0000000..e65216e
--- /dev/null
+++ b/man/initialize.FUN.Rd
@@ -0,0 +1,21 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/epi-functions.R
+\name{initialize.FUN}
+\alias{initialize.FUN}
+\title{Initialize ICM}
+\usage{
+initialize.FUN(param, init, control)
+}
+\arguments{
+\item{param}{ICM parameters.}
+
+\item{init}{Initial value parameters.}
+
+\item{control}{Control parameters.}
+}
+\value{
+Updated dat
+}
+\description{
+Function to initialize the icm
+}
diff --git a/man/merge.seiqhrf.icm.Rd b/man/merge.seiqhrf.icm.Rd
new file mode 100644
index 0000000..7bcc4b0
--- /dev/null
+++ b/man/merge.seiqhrf.icm.Rd
@@ -0,0 +1,19 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/sirPlus-functions.R
+\name{merge.seiqhrf.icm}
+\alias{merge.seiqhrf.icm}
+\title{Merge icm}
+\usage{
+\method{merge}{seiqhrf.icm}(x, y, ...)
+}
+\arguments{
+\item{x}{First to merge}
+
+\item{y}{Second to merge}
+}
+\value{
+Combined
+}
+\description{
+Function to merge icms
+}
diff --git a/man/progress.seiqhrf.icm.Rd b/man/progress.seiqhrf.icm.Rd
new file mode 100644
index 0000000..2ffb896
--- /dev/null
+++ b/man/progress.seiqhrf.icm.Rd
@@ -0,0 +1,19 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/sirPlus-functions.R
+\name{progress.seiqhrf.icm}
+\alias{progress.seiqhrf.icm}
+\title{Progress icm}
+\usage{
+progress.seiqhrf.icm(dat, at)
+}
+\arguments{
+\item{dat}{Object containing all data}
+
+\item{at}{?}
+}
+\value{
+progress
+}
+\description{
+Function to get progress of icms
+}
diff --git a/man/saveout.seiqhrf.icm.Rd b/man/saveout.seiqhrf.icm.Rd
new file mode 100644
index 0000000..31c337b
--- /dev/null
+++ b/man/saveout.seiqhrf.icm.Rd
@@ -0,0 +1,21 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/sirPlus-functions.R
+\name{saveout.seiqhrf.icm}
+\alias{saveout.seiqhrf.icm}
+\title{Save icm}
+\usage{
+saveout.seiqhrf.icm(dat, s, out = NULL)
+}
+\arguments{
+\item{dat}{Object containing all data}
+
+\item{out}{?}
+
+\item{at}{?}
+}
+\value{
+out to save
+}
+\description{
+Function to save icm
+}
diff --git a/man/set.control.Rd b/man/set.control.Rd
new file mode 100644
index 0000000..45c87fc
--- /dev/null
+++ b/man/set.control.Rd
@@ -0,0 +1,54 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/epi-functions.R
+\name{set.control}
+\alias{set.control}
+\title{Set control params}
+\usage{
+set.control(
+ 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
+)
+}
+\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{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.}
+}
+\value{
+control.icm object
+}
+\description{
+Sets the controls for stochastic individual contact models simulated with
+icm.
+}
diff --git a/man/set.init.Rd b/man/set.init.Rd
new file mode 100644
index 0000000..b649bbc
--- /dev/null
+++ b/man/set.init.Rd
@@ -0,0 +1,47 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/epi-functions.R
+\name{set.init}
+\alias{set.init}
+\title{Set init params}
+\usage{
+set.init(
+ s.num = 9997,
+ e.num = 0,
+ i.num = 3,
+ q.num = 0,
+ h.num = 0,
+ r.num = 0,
+ f.num = 0
+)
+}
+\arguments{
+\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.}
+}
+\value{
+init.icm object
+}
+\description{
+Sets the initial conditions for stochastic individual contact models
+simulated with icm.
+}
diff --git a/man/set.param.Rd b/man/set.param.Rd
new file mode 100644
index 0000000..33846aa
--- /dev/null
+++ b/man/set.param.Rd
@@ -0,0 +1,175 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/epi-functions.R
+\name{set.param}
+\alias{set.param}
+\title{Set epidemic params}
+\usage{
+set.param(
+ 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 = 1/20,
+ 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{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{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{
+param.icm object
+}
+\description{
+Sets the epidemic parameters for stochastic individual contact models
+simulated with icm.
+}
diff --git a/man/sirPlus.Rd b/man/sirPlus.Rd
new file mode 100644
index 0000000..0c18f8b
--- /dev/null
+++ b/man/sirPlus.Rd
@@ -0,0 +1,10 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/sirPlus-package.R
+\docType{package}
+\name{sirPlus}
+\alias{sirPlus}
+\title{sirPlus}
+\description{
+\pkg{sirPlus} is a package for modeling COVID-19 spread, focusing on rates
+of hospitalization, which should be useful for hospitals.
+}
diff --git a/vignettes/epimodels-SEIQHRF.Rmd b/vignettes/epimodels-SEIQHRF.Rmd
new file mode 100644
index 0000000..6c7ffad
--- /dev/null
+++ b/vignettes/epimodels-SEIQHRF.Rmd
@@ -0,0 +1,133 @@
+---
+title: "sirPLUS models"
+date: "Last updated: 23 March 2020"
+output:
+ BiocStyle::html_document:
+ toc: true
+ toc_float: true
+vignette: >
+ %\VignetteIndexEntry{sirPLUS models}
+ %\VignetteEngine{knitr::rmarkdown}
+ %\VignetteEncoding{UTF-8}
+---
+
+```{r, include = FALSE}
+knitr::opts_chunk$set(
+ collapse = TRUE,
+ comment = "#>"
+)
+```
+
+## Backgrouund on the SEIQHRF (SIR + extra compartments) model
+
+Based on [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/).
+
+
+*Note the lower case letters between compartment nodes represent model input
+parameters defined in the blog post*
+
+
+## Sett up environment and pull parameters
+
+```{r setup environment, include=FALSE}
+library(tidyverse)
+library(dplyr)
+library(magrittr)
+library(lubridate)
+library(stringr)
+library(tibble)
+library(broom)
+library(ggplot2)
+#remotes::install_github("rstudio/gt")
+library(gt)
+library(knitr)
+library(devtools)
+library(DiagrammeR)
+library(parallel)
+library(foreach)
+library(tictoc)
+suppressMessages(library(EpiModel))
+library(incidence)
+.libPaths('/mnt/mcfiles/rlyu/Software/R/3.6/Rlib')
+#devtools::install_github("reconhub/earlyR")
+library(earlyR)
+
+tic("Time to complete")
+```
+
+```{r data from Tim Churches}
+source_files <- c("_icm.mod.init.seiqhrf.R", "_icm.mod.status.seiqhrf.R",
+ "_icm.mod.vital.seiqhrf.R", "_icm.control.seiqhrf.R", "_icm.utils.seiqhrf.R",
+ "_icm.saveout.seiqhrf.R", "_icm.icm.seiqhrf.R")
+
+src_path <- paste0("./_posts/2020-03-18-modelling-the-effects-of-public-health-",
+ "interventions-on-covid-19-transmission-part-2/")
+
+gist_url <- "https://gist.github.com/timchurches/92073d0ea75cfbd387f91f7c6e624bd7"
+
+local_source <- FALSE
+
+for (source_file in source_files) {
+ if (local_source) {
+ source(paste(src_path, source_file, sep = ""))
+ } else {
+ Sys.sleep(1)
+ source_gist(gist_url, filename = source_file)
+ }
+}
+```
+
+## Define Functions
+
+```{r simulate baselines}
+# function to set-up and run the baseline simulations
+devtools::load_all(".")
+control <- set.control()
+param <- set.param()
+init <- set.init(s.num = 10000, e.num = 150, q.num = 10, h.num = 0)
+init
+
+sim <- icm.seiqhrf(param, init, control)
+sim_df <- as.data.frame(sim, out=out)
+simulate <- list(sim=sim, df=sim_df)
+
+```
+
+
+## Generate and inspect baseline simulations
+started at 3:12 pm
+```{r baseline sims}
+baseline_sim <- simulate(ncores = 4)
+times <- get_times(baseline_sim)
+
+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")
+```
+
+
+
+```{r viz prevalence}
+baseline_plot_df <- baseline_sim$df %>% # use only the prevalence columns
+select(time, s.num, e.num, i.num, q.num, h.num, r.num, f.num) %>%
+ # examine only the first 100 days since it is all over by
+# then using the default parameters
+filter(time <= 100) %>% pivot_longer(-c(time), names_to = "compartment",
+ values_to = "count")
+
+# define a standard set of colours to represent compartments
+compcols <- c(s.num = "yellow", e.num = "orange", i.num = "red",
+ q.num = "cyan", h.num = "magenta", r.num = "lightgreen",
+ f.num = "black")
+complabels <- c(s.num = "Susceptible", e.num = "Infected/asymptomatic",
+ i.num = "Infected/infectious", q.num = "Self-isolated", h.num = "Requires hospitalisation",
+ r.num = "Recovered", f.num = "Case fatality")
+
+baseline_plot_df %>% ggplot(aes(x = time, y = count, colour = compartment)) +
+ geom_line(size = 2, alpha = 0.7) + scale_colour_manual(values = compcols,
+ labels = complabels) + theme_dark() + labs(title = "Baseline simulation",
+ x = "Days since beginning of epidemic", y = "Prevalence (persons)")
+```
+## Experiment 1
+
+
--
GitLab