From 326ed0c37238ac283fc9ad3c83a4884996cd40b0 Mon Sep 17 00:00:00 2001
From: pqiao29 <pqiao@student.unimelb.edu.au>
Date: Mon, 6 Apr 2020 17:37:01 +1000
Subject: [PATCH] add from E to R: arec.rand, arec.rate, arec.dist.shape,
arec.dist.scale
---
R/FUN_recovery.R | 18 +++-
R/set_param.R | 135 +++++++++++++++++-------------
R/simulation.R | 67 +++++++++++++--
man/simulate_seiqhrf.Rd | 102 +++++++++++++++-------
tests/testthat/test-icm-seiqhrf.R | 2 +-
5 files changed, 229 insertions(+), 95 deletions(-)
diff --git a/R/FUN_recovery.R b/R/FUN_recovery.R
index 09c7d26..dfce50a 100644
--- a/R/FUN_recovery.R
+++ b/R/FUN_recovery.R
@@ -125,7 +125,23 @@ recovery.FUN <- function(dat, at, seed = NULL) {
dat$attr$recovTime[res[[2]]] <- at
status <- res[[3]]
dat$attr$status <- status
- # cat("finished recover \n")
+
+
+ res <- update_status(rate = dat$param$arec.rate,
+ rand = dat$control$arec.rand,
+ active,
+ status = dat$attr$status,
+ label = c("e"),
+ state = recovState,
+ at,
+ expTime = dat$attr$expTime,
+ prog.dist.scale = dat$param$arec.dist.scale,
+ prog.dist.shape = dat$param$arec.dist.shape)
+
+ nRecov <- nRecov + res[[1]]
+ dat$attr$recovTime[res[[2]]] <- at
+ status <- res[[3]]
+ dat$attr$status <- status
if (type %in% c("SEIQHRF")) {
# -------------------------------------- case fatality ----------------------------------------
diff --git a/R/set_param.R b/R/set_param.R
index 1731d6d..e1e79af 100644
--- a/R/set_param.R
+++ b/R/set_param.R
@@ -1,60 +1,70 @@
### internal function
set_param <- 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,
- infection.FUN = 'infection.FUN', # infection.seiqhrf.icm,
- recovery.FUN = 'recovery.FUN', # progress.seiqhrf.icm,
- departures.FUN = 'departures.FUN', # departures.seiqhrf.icm,
- arrivals.FUN = 'arrivals.FUN', # arrivals.seiqhrf.icm,
- get_prev.FUN = 'get_prev.FUN', # get_prev.seiqhrf.icm,
- # init.icm params
- s.num = 9997,
- e.num=0,
- i.num = 3,
- q.num=0,
- h.num=0,
- r.num = 0,
- f.num = 0,
- # param.icm params
- inf.prob.e = 0.02,
- act.rate.e = 10,
- inf.prob.i = 0.05,
- act.rate.i = 10,
- inf.prob.q = 0.02,
- act.rate.q = 2.5,
- quar.rate = 1/30,
- hosp.rate = 1/100,
- disch.rate = 1/15,
- prog.rate = 1/10,
- prog.dist.scale = 5,
- prog.dist.shape = 1.5,
- rec.rate = 0.071, # Based on assumption that disease course is 14 days
- 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"
+ nsteps = 366,
+ nsims = 8,
+ ncores = 4,
+ prog.rand = FALSE,
+ quar.rand = TRUE,
+ hosp.rand = TRUE,
+ disch.rand = TRUE,
+ rec.rand = TRUE,
+ arec.rand = TRUE,
+ fat.rand = TRUE,
+ infection.FUN = 'infection.FUN', # infection.seiqhrf.icm,
+ recovery.FUN = 'recovery.FUN', # progress.seiqhrf.icm,
+ departures.FUN = 'departures.FUN', # departures.seiqhrf.icm,
+ arrivals.FUN = 'arrivals.FUN', # arrivals.seiqhrf.icm,
+ get_prev.FUN = 'get_prev.FUN', # get_prev.seiqhrf.icm,
+ # init.icm params
+ s.num = 9997,
+ e.num=0,
+ i.num = 3,
+ q.num=0,
+ h.num=0,
+ r.num = 0,
+ f.num = 0,
+ # param.icm params
+ 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,
+ prog.rate = 1/10,
+ quar.rate = 1/30,
+ hosp.rate = 1/100,
+ disch.rate = 1/15,
+ rec.rate = 0.071,
+ arec.rate = 0.05,
+ prog.dist.scale = 5,
+ prog.dist.shape = 1.5,
+ quar.dist.scale = 1,
+ quar.dist.shape = 1,
+ hosp.dist.scale = 1,
+ hosp.dist.shape = 1,
+ disch.dist.scale = 1,
+ disch.dist.shape = 1,
+ rec.dist.scale = 35,
+ rec.dist.shape = 1.5,
+ arec.dist.scale = 35,
+ arec.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"
) {
control <- control.icm(type = type,
@@ -62,7 +72,11 @@ set_param <- function(type = "SEIQHRF",
nsims = nsims,
ncores = ncores,
prog.rand = prog.rand,
+ quar.rand = quar.rand,
+ hosp.ramd = hosp.rand,
+ disch.rand = disch.rand,
rec.rand = rec.rand,
+ arec.rand = arec.rand,
infection.FUN = infection.FUN,
recovery.FUN = recovery.FUN,
arrivals.FUN = arrivals.FUN,
@@ -83,15 +97,24 @@ set_param <- function(type = "SEIQHRF",
act.rate.i = act.rate.i,
inf.prob.q = inf.prob.q,
act.rate.q = act.rate.q,
+ prog.rate = prog.rate,
quar.rate = quar.rate,
hosp.rate = hosp.rate,
disch.rate = disch.rate,
- prog.rate = prog.rate,
+ rec.rate = rec.rate,
+ arec.rate = arec.rate,
prog.dist.scale = prog.dist.scale,
prog.dist.shape = prog.dist.shape,
- rec.rate = rec.rate,
+ quar.dist.scale = quar.dist.scale,
+ quar.dist.shape = quar.dist.shape,
+ hosp.dist.scale = hosp.dist.scale,
+ hosp.dist.shape = hosp.dist.shape,
+ disch.dist.scale = disch.dist.scale,
+ disch.dist.shape = disch.dist.shape,
rec.dist.scale = rec.dist.scale,
rec.dist.shape = rec.dist.shape,
+ arec.dist.scale = arec.dist.scale,
+ arec.dist.shape = arec.dist.shape,
fat.rate.base = fat.rate.base,
hosp.cap = hosp.cap,
fat.rate.overcap = fat.rate.overcap,
diff --git a/R/simulation.R b/R/simulation.R
index 29bcb52..630170b 100644
--- a/R/simulation.R
+++ b/R/simulation.R
@@ -15,9 +15,13 @@
#' distribution but it makes little difference and speed of computation
#' matters), with parameters prog.dist.scale and prog.dist.shape
#' @param rec.rand Method for recovery transition from I, Q or H to R. If TRUE,
-#' random binomial draws at prog.rate, if FALSE, random draws from a
+#' random binomial draws at rec.rate, if FALSE, random draws from a
#' random draws from a Weibull distribution, with parameters
#' rec.dist.scale and rec.dist.shape.
+#' @param arec.rand Method for recovery transition from E to R. If TRUE,
+#' random binomial draws at arec.rate, if FALSE, random draws from a
+#' random draws from a Weibull distribution, with parameters
+#' arec.dist.scale and raec.dist.shape.
#' @param fat.rand Method for case fatality transition from H to F. If TRUE,
#' random binomial draws at fat.rate.base, if FALSE, random sample with
#' a sample fraction also given by fat.rate.base. However, if the
@@ -117,6 +121,11 @@
#' a low rate reflecting low community awareness or compliance with
#' self-isolation requirements or practices, but this can be tweaked when
#' exploring scenarios.
+#' @param quar.dist.scale Scale parameter for Weibull distribution for
+#' recovery, see quar.rand for details.
+#' @param quar.dist.shape Shape parameter for Weibull distribution for
+#' recovery, see quar.rand for details. Read up on the Weibull distribution
+#' before changing the default.
#' @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
@@ -124,8 +133,18 @@
#' 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 hosp.dist.scale Scale parameter for Weibull distribution for
+#' recovery, see hosp.rand for details.
+#' @param hosp.dist.shape Shape parameter for Weibull distribution for
+#' recovery, see hosp.rand for details. Read up on the Weibull distribution
+#' before changing the default.
#' @param disch.rate Rate per day at which people needing
#' hospitalisation recover.
+#' @param disch.dist.scale Scale parameter for Weibull distribution for
+#' recovery, see disch.rand for details.
+#' @param disch.dist.shape Shape parameter for Weibull distribution for
+#' recovery, see disch.rand for details. Read up on the Weibull distribution
+#' before changing the default.
#' @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.
@@ -142,6 +161,13 @@
#' @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 arec.rate Rate per day at which people who are exposed but asymptotic
+#' (E compartment) recover, thus entering the R compartment.
+#' See arec.rand above for more details.
+#' @param arec.dist.scale Scale parameter for Weibull distribution for
+#' recovery, see raec.rand for details.
+#' @param arec.dist.shape Shape parameter for Weibull distribution for
+#' recovery, see arec.rand for details.
#' @param fat.rate.base Baseline mortality rate per day for people
#' needing hospitalisation (deaths due to the virus). See fat.rand for more
#' details.
@@ -193,11 +219,12 @@ simulate_seiqhrf <- function(type = "SEIQHRF",
nsims = 8,
ncores = 4,
prog.rand = FALSE,
- rec.rand = FALSE,
- fat.rand = TRUE,
- quar.rand = FALSE,
- hosp.rand = FALSE,
+ quar.rand = TRUE,
+ hosp.rand = TRUE,
disch.rand = TRUE,
+ rec.rand = TRUE,
+ arec.rand = TRUE,
+ fat.rand = TRUE,
infection.FUN = 'infection.FUN', # infection.seiqhrf.icm,
recovery.FUN = 'recovery.FUN', # progress.seiqhrf.icm,
departures.FUN = 'departures.FUN', # departures.seiqhrf.icm,
@@ -218,15 +245,24 @@ simulate_seiqhrf <- function(type = "SEIQHRF",
act.rate.i = 10,
inf.prob.q = 0.02,
act.rate.q = 2.5,
+ prog.rate = 1/10,
quar.rate = 1/30,
hosp.rate = 1/100,
disch.rate = 1/15,
- prog.rate = 1/10,
+ rec.rate = 0.071,
+ arec.rate = 0.05,
prog.dist.scale = 5,
prog.dist.shape = 1.5,
- rec.rate = 0.071, # Based on assumption that disease course is 14 days
+ quar.dist.scale = 1,
+ quar.dist.shape = 1,
+ hosp.dist.scale = 1,
+ hosp.dist.shape = 1,
+ disch.dist.scale = 1,
+ disch.dist.shape = 1,
rec.dist.scale = 35,
rec.dist.shape = 1.5,
+ arec.dist.scale = 35,
+ arec.dist.shape = 1.5,
fat.rate.base = 1/50,
hosp.cap = 40,
fat.rate.overcap = 1/25,
@@ -250,7 +286,11 @@ simulate_seiqhrf <- function(type = "SEIQHRF",
nsims = nsims,
ncores = ncores,
prog.rand = prog.rand,
+ quar.rand = quar.rand,
+ hosp.ramd = hosp.rand,
+ disch.rand = disch.rand,
rec.rand = rec.rand,
+ arec.rand = arec.rand,
infection.FUN = infection.FUN,
recovery.FUN = recovery.FUN,
arrivals.FUN = arrivals.FUN,
@@ -271,15 +311,24 @@ simulate_seiqhrf <- function(type = "SEIQHRF",
act.rate.i = act.rate.i,
inf.prob.q = inf.prob.q,
act.rate.q = act.rate.q,
+ prog.rate = prog.rate,
quar.rate = quar.rate,
hosp.rate = hosp.rate,
disch.rate = disch.rate,
- prog.rate = prog.rate,
+ rec.rate = rec.rate,
+ arec.rate = arec.rate,
prog.dist.scale = prog.dist.scale,
prog.dist.shape = prog.dist.shape,
- rec.rate = rec.rate,
+ quar.dist.scale = quar.dist.scale,
+ quar.dist.shape = quar.dist.shape,
+ hosp.dist.scale = hosp.dist.scale,
+ hosp.dist.shape = hosp.dist.shape,
+ disch.dist.scale = disch.dist.scale,
+ disch.dist.shape = disch.dist.shape,
rec.dist.scale = rec.dist.scale,
rec.dist.shape = rec.dist.shape,
+ arec.dist.scale = arec.dist.scale,
+ arec.dist.shape = arec.dist.shape,
fat.rate.base = fat.rate.base,
hosp.cap = hosp.cap,
fat.rate.overcap = fat.rate.overcap,
diff --git a/man/simulate_seiqhrf.Rd b/man/simulate_seiqhrf.Rd
index 7d246a4..e8f014d 100644
--- a/man/simulate_seiqhrf.Rd
+++ b/man/simulate_seiqhrf.Rd
@@ -10,11 +10,12 @@ simulate_seiqhrf(
nsims = 8,
ncores = 4,
prog.rand = FALSE,
- rec.rand = FALSE,
- fat.rand = TRUE,
- quar.rand = FALSE,
- hosp.rand = FALSE,
+ quar.rand = TRUE,
+ hosp.rand = TRUE,
disch.rand = TRUE,
+ rec.rand = TRUE,
+ arec.rand = TRUE,
+ fat.rand = TRUE,
infection.FUN = "infection.FUN",
recovery.FUN = "recovery.FUN",
departures.FUN = "departures.FUN",
@@ -33,15 +34,24 @@ simulate_seiqhrf(
act.rate.i = 10,
inf.prob.q = 0.02,
act.rate.q = 2.5,
+ prog.rate = 1/10,
quar.rate = 1/30,
hosp.rate = 1/100,
disch.rate = 1/15,
- prog.rate = 1/10,
+ rec.rate = 0.071,
+ arec.rate = 0.05,
prog.dist.scale = 5,
prog.dist.shape = 1.5,
- rec.rate = 0.071,
+ quar.dist.scale = 1,
+ quar.dist.shape = 1,
+ hosp.dist.scale = 1,
+ hosp.dist.shape = 1,
+ disch.dist.scale = 1,
+ disch.dist.shape = 1,
rec.dist.scale = 35,
rec.dist.shape = 1.5,
+ arec.dist.scale = 35,
+ arec.dist.shape = 1.5,
fat.rate.base = 1/50,
hosp.cap = 40,
fat.rate.overcap = 1/25,
@@ -78,11 +88,31 @@ Weibull distribution (yes, I know it should be a discrete Weibull
distribution but it makes little difference and speed of computation
matters), with parameters prog.dist.scale and prog.dist.shape}
+\item{quar.rand}{Method for self-isolation transition from I to Q. If TRUE,
+random binomial draws at quar.rate, if FALSE, random sample with a
+sample fraction also given by `quar.rate.}
+
+\item{hosp.rand}{Method for transition from I or Q to H -- that is, from
+infectious or from self-isolated to requiring hospitalisation. If
+TRUE, random binomial draws at hosp.rate, if FALSE, random sample
+with a sample fraction also given by `hosp.rate.}
+
+\item{disch.rand}{Method for transition from H to R -- that is, from
+requiring hospitalisation to recovered. If TRUE, random binomial
+draws at disch.rate, if FALSE, random sample with a sample fraction
+also given by disch.rate. Note that the only way out of the H
+compartment is recovery or death.}
+
\item{rec.rand}{Method for recovery transition from I, Q or H to R. If TRUE,
-random binomial draws at prog.rate, if FALSE, random draws from a
+random binomial draws at rec.rate, if FALSE, random draws from a
random draws from a Weibull distribution, with parameters
rec.dist.scale and rec.dist.shape.}
+\item{arec.rand}{Method for recovery transition from E to R. If TRUE,
+random binomial draws at arec.rate, if FALSE, random draws from a
+random draws from a Weibull distribution, with parameters
+arec.dist.scale and raec.dist.shape.}
+
\item{fat.rand}{Method for case fatality transition from H to F. If TRUE,
random binomial draws at fat.rate.base, if FALSE, random sample with
a sample fraction also given by fat.rate.base. However, if the
@@ -99,21 +129,6 @@ as a linear function of the number of days the individual has been
in the H compartment. Use of the co-efficient better approximates
the trapezoid survival time distribution typical of ICU patients.}
-\item{quar.rand}{Method for self-isolation transition from I to Q. If TRUE,
-random binomial draws at quar.rate, if FALSE, random sample with a
-sample fraction also given by `quar.rate.}
-
-\item{hosp.rand}{Method for transition from I or Q to H -- that is, from
-infectious or from self-isolated to requiring hospitalisation. If
-TRUE, random binomial draws at hosp.rate, if FALSE, random sample
-with a sample fraction also given by `hosp.rate.}
-
-\item{disch.rand}{Method for transition from H to R -- that is, from
-requiring hospitalisation to recovered. If TRUE, random binomial
-draws at disch.rate, if FALSE, random sample with a sample fraction
-also given by disch.rate. Note that the only way out of the H
-compartment is recovery or death.}
-
\item{infection.FUN}{No, being infected with SARS-CoV2 is not fun. Rather
this is the the name of the function to implement infection
processes. Use the default.}
@@ -197,6 +212,10 @@ 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{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{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
@@ -216,9 +235,13 @@ 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{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{arec.rate}{Rate per day at which people who are exposed but asymptotic
+(E compartment) recover, thus entering the R compartment.
+See arec.rand above for more details.}
\item{prog.dist.scale}{Scale parameter for Weibull distribution
for progression, see prog.rand for details.}
@@ -227,9 +250,26 @@ for progression, see prog.rand for details.}
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{quar.dist.scale}{Scale parameter for Weibull distribution for
+recovery, see quar.rand for details.}
+
+\item{quar.dist.shape}{Shape parameter for Weibull distribution for
+recovery, see quar.rand for details. Read up on the Weibull distribution
+before changing the default.}
+
+\item{hosp.dist.scale}{Scale parameter for Weibull distribution for
+recovery, see hosp.rand for details.}
+
+\item{hosp.dist.shape}{Shape parameter for Weibull distribution for
+recovery, see hosp.rand for details. Read up on the Weibull distribution
+before changing the default.}
+
+\item{disch.dist.scale}{Scale parameter for Weibull distribution for
+recovery, see disch.rand for details.}
+
+\item{disch.dist.shape}{Shape parameter for Weibull distribution for
+recovery, see disch.rand for details. Read up on the Weibull distribution
+before changing the default.}
\item{rec.dist.scale}{Scale parameter for Weibull distribution for
recovery, see rec.rand for details.}
@@ -238,6 +278,12 @@ recovery, see rec.rand for details.}
recovery, see rec.rand for details. Read up on the Weibull distribution
before changing the default.}
+\item{arec.dist.scale}{Scale parameter for Weibull distribution for
+recovery, see raec.rand for details.}
+
+\item{arec.dist.shape}{Shape parameter for Weibull distribution for
+recovery, see arec.rand for details.}
+
\item{fat.rate.base}{Baseline mortality rate per day for people
needing hospitalisation (deaths due to the virus). See fat.rand for more
details.}
diff --git a/tests/testthat/test-icm-seiqhrf.R b/tests/testthat/test-icm-seiqhrf.R
index b355d8e..e9edfd6 100644
--- a/tests/testthat/test-icm-seiqhrf.R
+++ b/tests/testthat/test-icm-seiqhrf.R
@@ -1,6 +1,6 @@
test_that("Identical output as Churches' original function: icm_seiqhrf", {
- full_params <- set_param(s.num = 1000, nsteps = 10)
+ full_params <- set_param(s.num = 1000, nsteps = 10, arec.rate = 0)
param <- full_params$param
init <- full_params$init
--
GitLab