diff --git a/DESCRIPTION b/DESCRIPTION
index cb9f38c681feed5761007decbcc99256a7901d0d..d0ce1efca49f092268003da33ecf94b250c464a5 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -39,7 +39,8 @@ Suggests:
magick,
rmarkdown,
spelling,
- testthat
+ testthat,
+ devtools
biocViews: Epidemiology, Software
RoxygenNote: 7.1.0
VignetteBuilder: knitr
diff --git a/R/recovery.FUN.R b/R/recovery.FUN.R
index eadc963886a13d83c25b0c90ccc913e76249834f..37af7406c7364ecce640893c8c92176c95ef30c5 100644
--- a/R/recovery.FUN.R
+++ b/R/recovery.FUN.R
@@ -4,7 +4,7 @@ update_status <- function(rate, rand, active, status, label, state, at, prog, ex
smp_sz <- 0
at_idx <- NULL
- idsElig <- which(active == 1 & status == label)
+ idsElig <- which(active == 1 & status %in% label)
nElig <- length(idsElig)
if (nElig > 0) {
@@ -13,10 +13,9 @@ update_status <- function(rate, rand, active, status, label, state, at, prog, ex
if (rand) {
vecProg <- which(stats::rbinom(nElig, 1, rate) == 1)
if (length(vecProg) > 0) {
- idsProg <- idsElig[vecProg]
- smp_sz <- length(idsProg)
- status[idsProg] <- state
- at_idx <- idsProg
+ at_idx <- idsElig[vecProg]
+ smp_sz <- length(at_idx)
+ status[at_idx] <- state
}
}else{
do_sample <- TRUE
@@ -25,6 +24,7 @@ update_status <- function(rate, rand, active, status, label, state, at, prog, ex
smp_prob <- NULL
}else{
vecTimeSinceExp <- at - expTime[idsElig]
+ vecTimeSinceExp[is.na(vecTimeSinceExp)] <- 0
gammaRatesElig <- stats::pweibull(vecTimeSinceExp, prog.dist.shape, scale=prog.dist.scale)
smp_sz <- round(sum(gammaRatesElig[gElig == 1], na.rm=TRUE))
smp_prob <- gammaRatesElig[gElig == 1]
@@ -32,9 +32,8 @@ update_status <- function(rate, rand, active, status, label, state, at, prog, ex
}
if(do_sample){
- ids2bProg <- EpiModel::ssample(idsElig[gElig == 1], smp_sz, prob = smp_prob)
- status[ids2bProg] <- state
- at_idx <- ids2bProg
+ at_idx <- EpiModel::ssample(idsElig[gElig == 1], smp_sz, prob = smp_prob)
+ status[at_idx] <- state
}
@@ -63,13 +62,15 @@ update_status <- function(rate, rand, active, status, label, state, at, prog, ex
#' @export
recovery.FUN <- function(dat, at, seed = NULL) {
+ #cat("at: ", at, "\n")
+
if(!is.null(seed)) set.seed(seed)
- # Conditions --------------------------------------------------------------
+ # Conditions --------------------------------------------------------------------------
if (!(dat$control$type %in% c("SIR", "SIS", "SEIR", "SEIQHR", "SEIQHRF"))) {
return(dat)
}
- # Variables ---------------------------------------------------------------
+ # Variables ---------------------------------------------------------------------------
active <- dat$attr$active
status <- dat$attr$status
@@ -81,10 +82,10 @@ recovery.FUN <- function(dat, at, seed = NULL) {
hospState <- "h"
fatState <- "f"
- # --- progress from exposed to infectious ----
+ # ------------------------- progress from exposed to infectious --------------------------
res <- update_status(rate = dat$param$prog.rate,
rand = dat$control$prog.rand,
- active, status,
+ active, dat$attr$status,
label = "e",
state = "i",
at, prog = TRUE,
@@ -97,6 +98,8 @@ recovery.FUN <- function(dat, at, seed = NULL) {
status <- res[[3]]
dat$attr$status <- status
+ #cat("finished prog \n")
+
if (type %in% c("SEIQHR", "SEIQHRF")) {
# ----- quarantine -------
@@ -114,242 +117,110 @@ recovery.FUN <- function(dat, at, seed = NULL) {
status <- res[[3]]
dat$attr$status <- status
+ #cat("finished quarantine \n")
+
# ----- need to be hospitalised -------
- hosp.rand <- dat$control$hosp.rand
- hosp.rate <- dat$param$hosp.rate
+ rate <- rep(dat$param$hosp.rate,
+ sum(dat$attr$active == 1 & (dat$attr$status == "i" | dat$attr$status == "q")))
- nHosp <- nHospG2 <- 0
- idsElig <- which(active == 1 & (status == "i" | status == "q"))
- nElig <- length(idsElig)
- idsHosp <- numeric(0)
+ res <- update_status(rate,
+ rand = dat$control$hosp.rand,
+ active = dat$attr$active,
+ status = dat$attr$status,
+ label = c("i", "q"),
+ state = "h", at, prog = FALSE)
- if (nElig > 0) {
-
- gElig <- group[idsElig]
- rates <- hosp.rate
- ratesElig <- rates[gElig]
-
- if (hosp.rand) {
- vecHosp <- which(stats::rbinom(nElig, 1, ratesElig) == 1)
- if (length(vecHosp) > 0) {
- idsHosp <- idsElig[vecHosp]
- nHosp <- sum(group[idsHosp] == 1)
- status[idsHosp] <- hospState
- }
- } else {
- nHosp <- min(round(sum(ratesElig[gElig == 1])), sum(gElig == 1))
- idsHosp <- ssample(idsElig[gElig == 1], nHosp)
- status[idsHosp] <- hospState
- }
- }
+ nHosp <- res[[1]]
+ dat$attr$hospTime[res[[2]]] <- at
+ status <- res[[3]]
dat$attr$status <- status
- dat$attr$hospTime[idsHosp] <- at
- # ----- discharge from need to be hospitalised -------
- disch.rand <- dat$control$disch.rand
+ # cat("finished hospitalised \n")
+
+ # ------------------------- discharge from need to be hospitalised ---------------------------
+ recovState <- ifelse(type %in% c("SIR", "SEIR", "SEIQHR", "SEIQHRF"), "r", "s")
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)
- }
- }
- }
+ dcrate <- ifelse(length(disch.rate) > 1, disch.rate[at], disch.rate)
+
+ res <- update_status(rate = dcrate,
+ rand = dat$control$disch.rand,
+ active = dat$attr$active,
+ status = dat$attr$status,
+ label = "h",
+ state = recovState, at, prog = FALSE)
+
+ nDisch <- res[[1]]
+ dat$attr$dischTime[res[[2]]] <- at
+ status <- res[[3]]
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)
- }
- }
- }
+ #cat("finished discharge \n")
}
+
+ # --------------------------------------- recover -------------------------------------------------
+ res <- update_status(rate = dat$param$rec.rate,
+ rand = dat$control$rec.rand,
+ active,
+ status = dat$attr$status,
+ label = c("i", "q", "h"),
+ state = recovState,
+ at, prog = TRUE,
+ expTime = dat$attr$expTime,
+ prog.dist.scale = dat$param$rec.dist.scale,
+ prog.dist.shape = dat$param$rec.dist.shape)
+
+ nRecov <- res[[1]]
+ dat$attr$recovTime[res[[2]]] <- at
+ status <- res[[3]]
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
+ # cat("finished recover \n")
+
+ if (type %in% c("SEIQHRF")) {
+ # -------------------------------------- case fatality ----------------------------------------
+ ## obtain rates
idsElig <- which(active == 1 & status =="h")
- nElig <- length(idsElig)
-
- if (nElig > 0) {
+ if(length(idsElig) > 0){
+ hosp.cap <- dat$param$hosp.cap
+ rates <- dat$param$fat.rate.base
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)
+ blended.rate <- ((hosp.cap * rates) +
+ ((h.num.yesterday - hosp.cap) * dat$param$fat.rate.overcap)) /
+ h.num.yesterday
+ rates <- blended.rate
}
}
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
- }
- }
+ ratesElig <- ratesElig + timeInHospElig*dat$param$fat.tcoeff*ratesElig
}
+
+ res <- update_status(rate = ratesElig,
+ rand = dat$control$fat.rand,
+ active = dat$attr$active,
+ status = dat$attr$status,
+ label = "h",
+ state = "f", at, prog = FALSE)
+ nFat <- res[[1]]
+ if(!is.null(res[[2]])) dat$attr$fatTime[res[[2]]] <- at
+ status <- res[[3]]
dat$attr$status <- status
+
+ #cat("finished fatality \n")
}
-
+
# 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) {
@@ -361,7 +232,7 @@ recovery.FUN <- function(dat, at, seed = NULL) {
dat$epi[[outName_c[1]]] <- c(0, nQuar)
dat$epi[[outName_d[1]]] <- c(0, nHosp)
}
- if (fatEnable & type %in% c("SEIQHRF")) {
+ if (type %in% c("SEIQHRF")) {
dat$epi[[outName_e[1]]] <- c(0, nFat)
}
} else {
@@ -373,37 +244,10 @@ recovery.FUN <- function(dat, at, seed = NULL) {
dat$epi[[outName_c[1]]][at] <- nQuar
dat$epi[[outName_d[1]]][at] <- nHosp
}
- if (fatEnable & type %in% c("SEIQHRF")) {
+ if (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)
}