Skip to content
Snippets Groups Projects
Commit 4323e407 authored by pqiao29's avatar pqiao29
Browse files

finish clean up recovery.FUN

parent 6bc64cf5
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
DESCRIPTION
+ 2
− 1
View file @ 4323e407
......@@ -39,7 +39,8 @@ Suggests:
magick,
rmarkdown,
spelling,
testthat
testthat,
devtools
biocViews: Epidemiology, Software
RoxygenNote: 7.1.0
VignetteBuilder: knitr
......
R/recovery.FUN.R
+ 90
− 246
View file @ 4323e407
......@@ -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)
}
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment