1 |
#' Progress icm |
|
2 |
#' |
|
3 |
#' Function to get progress of icms |
|
4 |
#' |
|
5 |
#' @param dat Object containing all data |
|
6 |
#' @param at time point |
|
7 |
#' @param seed random seed for checking consistency |
|
8 |
#' |
|
9 |
#' @return progress |
|
10 |
#' @importFrom stats pweibull |
|
11 |
#' @importFrom stats rbinom |
|
12 |
#' @importFrom EpiModel ssample |
|
13 |
#' @importFrom stats rgeom |
|
14 |
#' @importFrom stats sd |
|
15 |
#' @export |
|
16 |
recovery.FUN <- function(dat, at, seed = NULL) { |
|
17 |
|
|
18 |
#cat("at: ", at, "\n") |
|
19 |
|
|
20 | ! |
if(!is.null(seed)) set.seed(seed) |
21 |
# Conditions -------------------------------------------------------------------------- |
|
22 | 2x |
if (!(dat$control$type %in% c("SIR", "SIS", "SEIR", "SEIQHR", "SEIQHRF"))) { |
23 | ! |
return(dat) |
24 |
} |
|
25 |
|
|
26 |
# Variables --------------------------------------------------------------------------- |
|
27 | 2x |
active <- dat$attr$active |
28 | 2x |
status <- dat$attr$status |
29 | 2x |
groups <- dat$param$groups |
30 | 2x |
group <- dat$attr$group |
31 | 2x |
type <- dat$control$type |
32 | 2x |
expTime <- dat$attr$expTime |
33 |
|
|
34 |
# ------------------------- progress from exposed to infectious -------------------------- |
|
35 | 2x |
res <- update_status(rate = dat$param$prog.rate, |
36 | 2x |
rand = dat$control$prog.rand, |
37 | 2x |
active, dat$attr$status, |
38 | 2x |
label = "e", |
39 | 2x |
state = "i", |
40 | 2x |
at, expTime, |
41 | 2x |
prog.dist.scale = dat$param$prog.dist.scale, |
42 | 2x |
prog.dist.shape = dat$param$prog.dist.shape) |
43 |
|
|
44 | 2x |
nProg <- res[[1]] |
45 | ! |
if(!is.null(res[[2]])) dat$attr$infTime[res[[2]]] <- at |
46 | 2x |
status <- res[[3]] |
47 | 2x |
dat$attr$status <- status |
48 |
|
|
49 |
#cat("finished prog \n") |
|
50 |
|
|
51 | 2x |
if (type %in% c("SEIQHR", "SEIQHRF")) { |
52 |
|
|
53 |
# ----- quarantine ------- |
|
54 | 2x |
quar.rate <- dat$param$quar.rate |
55 | 2x |
rate <- ifelse(length(quar.rate) > 1, quar.rate[at], quar.rate) |
56 |
|
|
57 | 2x |
res <- update_status(rate, |
58 | 2x |
rand = dat$control$quar.rand, |
59 | 2x |
active = dat$attr$active, |
60 | 2x |
status = dat$attr$status, |
61 | 2x |
label = "i", |
62 | 2x |
state = "q", |
63 | 2x |
at, expTime, |
64 | 2x |
prog.dist.scale = dat$param$quar.dist.scale, |
65 | 2x |
prog.dist.shape = dat$param$quar.dist.shape) |
66 | 2x |
nQuar <- res[[1]] |
67 | ! |
if(!is.null(res[[2]])) dat$attr$quarTime[res[[2]]] <- at |
68 | 2x |
status <- res[[3]] |
69 | 2x |
dat$attr$status <- status |
70 |
|
|
71 |
# ----- need to be hospitalised ------- |
|
72 | 2x |
rate <- rep(dat$param$hosp.rate, |
73 | 2x |
sum(dat$attr$active == 1 & (dat$attr$status == "i" | dat$attr$status == "q"))) |
74 |
|
|
75 | 2x |
res <- update_status(rate, |
76 | 2x |
rand = dat$control$hosp.rand, |
77 | 2x |
active = dat$attr$active, |
78 | 2x |
status = dat$attr$status, |
79 | 2x |
label = c("i", "q"), |
80 | 2x |
state = "h", |
81 | 2x |
at, expTime, |
82 | 2x |
prog.dist.scale = dat$param$hosp.dist.scale, |
83 | 2x |
prog.dist.shape = dat$param$hosp.dist.shape) |
84 |
|
|
85 | 2x |
nHosp <- res[[1]] |
86 | 2x |
dat$attr$hospTime[res[[2]]] <- at |
87 | 2x |
status <- res[[3]] |
88 | 2x |
dat$attr$status <- status |
89 |
|
|
90 |
# ------------------------- discharge from need to be hospitalised --------------------------- |
|
91 | 2x |
recovState <- ifelse(type %in% c("SIR", "SEIR", "SEIQHR", "SEIQHRF"), "r", "s") |
92 | 2x |
disch.rate <- dat$param$disch.rate |
93 | 2x |
dcrate <- ifelse(length(disch.rate) > 1, disch.rate[at], disch.rate) |
94 |
|
|
95 | 2x |
res <- update_status(rate = dcrate, |
96 | 2x |
rand = dat$control$disch.rand, |
97 | 2x |
active = dat$attr$active, |
98 | 2x |
status = dat$attr$status, |
99 | 2x |
label = "h", |
100 | 2x |
state = recovState, |
101 | 2x |
at, expTime, |
102 | 2x |
prog.dist.scale = dat$param$disch.dist.scale, |
103 | 2x |
prog.dist.shape = dat$param$disch.dist.shape) |
104 |
|
|
105 | 2x |
nDisch <- res[[1]] |
106 | 2x |
dat$attr$dischTime[res[[2]]] <- at |
107 | 2x |
status <- res[[3]] |
108 | 2x |
dat$attr$status <- status |
109 |
#cat("finished discharge \n") |
|
110 |
} |
|
111 |
|
|
112 |
# --------------------------------------- recover ------------------------------------------------- |
|
113 | 2x |
res <- update_status(rate = dat$param$rec.rate, |
114 | 2x |
rand = dat$control$rec.rand, |
115 | 2x |
active, |
116 | 2x |
status = dat$attr$status, |
117 | 2x |
label = c("i", "q"), |
118 | 2x |
state = recovState, |
119 | 2x |
at, |
120 | 2x |
expTime = dat$attr$expTime, |
121 | 2x |
prog.dist.scale = dat$param$rec.dist.scale, |
122 | 2x |
prog.dist.shape = dat$param$rec.dist.shape) |
123 |
|
|
124 | 2x |
nRecov <- res[[1]] |
125 | 2x |
dat$attr$recovTime[res[[2]]] <- at |
126 | 2x |
status <- res[[3]] |
127 | 2x |
dat$attr$status <- status |
128 |
|
|
129 |
|
|
130 | 2x |
res <- update_status(rate = dat$param$arec.rate, |
131 | 2x |
rand = dat$control$arec.rand, |
132 | 2x |
active, |
133 | 2x |
status = dat$attr$status, |
134 | 2x |
label = c("e"), |
135 | 2x |
state = recovState, |
136 | 2x |
at, |
137 | 2x |
expTime = dat$attr$expTime, |
138 | 2x |
prog.dist.scale = dat$param$arec.dist.scale, |
139 | 2x |
prog.dist.shape = dat$param$arec.dist.shape) |
140 |
|
|
141 | 2x |
nRecov <- nRecov + res[[1]] |
142 | 2x |
dat$attr$recovTime[res[[2]]] <- at |
143 | 2x |
status <- res[[3]] |
144 | 2x |
dat$attr$status <- status |
145 |
|
|
146 | 2x |
if (type %in% c("SEIQHRF")) { |
147 |
# -------------------------------------- case fatality ---------------------------------------- |
|
148 |
## obtain rates |
|
149 | 2x |
idsElig <- which(active == 1 & status =="h") |
150 | 2x |
if(length(idsElig) > 0){ |
151 | ! |
hosp.cap <- dat$param$hosp.cap |
152 | ! |
rates <- dat$param$fat.rate.base |
153 | ! |
gElig <- group[idsElig] |
154 | ! |
timeInHospElig <- at - dat$attr$hospTime[idsElig] |
155 | ! |
h.num.yesterday <- 0 |
156 | ! |
if (!is.null(dat$epi$h.num[at - 1])) { |
157 | ! |
h.num.yesterday <- dat$epi$h.num[at - 1] |
158 | ! |
if (h.num.yesterday > hosp.cap) { |
159 | ! |
blended.rate <- ((hosp.cap * rates) + |
160 | ! |
((h.num.yesterday - hosp.cap) * dat$param$fat.rate.overcap)) / |
161 | ! |
h.num.yesterday |
162 | ! |
rates <- blended.rate |
163 |
} |
|
164 |
} |
|
165 | ! |
ratesElig <- rates[gElig] |
166 | ! |
ratesElig <- ratesElig + timeInHospElig*dat$param$fat.tcoeff*ratesElig |
167 |
} |
|
168 |
|
|
169 | 2x |
res <- update_status(rate = ratesElig, |
170 | 2x |
rand = dat$control$fat.rand, |
171 | 2x |
active = dat$attr$active, |
172 | 2x |
status = dat$attr$status, |
173 | 2x |
label = "h", |
174 | 2x |
state = "f", at) |
175 | 2x |
nFat <- res[[1]] |
176 | ! |
if(!is.null(res[[2]])) dat$attr$fatTime[res[[2]]] <- at |
177 | 2x |
status <- res[[3]] |
178 | 2x |
dat$attr$status <- status |
179 |
|
|
180 |
#cat("finished fatality \n") |
|
181 |
} |
|
182 |
|
|
183 |
# Output ------------------------------------------------------------------ |
|
184 | 2x |
outName_a <- ifelse(type %in% c("SIR", "SEIR"), "ir.flow", "is.flow") |
185 | 2x |
if (type %in% c("SEIR", "SEIQHR", "SEIQHRF")) { |
186 | 2x |
outName_b <- "ei.flow" |
187 |
} |
|
188 | 2x |
if (type %in% c("SEIQHR", "SEIQHRF")) { |
189 | 2x |
outName_c <- "iq.flow" |
190 | 2x |
outName_d <- "iq2h.flow" |
191 |
} |
|
192 | 2x |
if (type %in% c("SEIQHRF")) { |
193 | 2x |
outName_e <- "hf.flow" |
194 |
} |
|
195 |
## Summary statistics |
|
196 | 2x |
if (at == 2) { |
197 | 2x |
dat$epi[[outName_a[1]]] <- c(0, nRecov) |
198 | 2x |
if (type %in% c("SEIR", "SEIQHR")) { |
199 | ! |
dat$epi[[outName_b[1]]] <- c(0, nProg) |
200 |
} |
|
201 | 2x |
if (type %in% c("SEIQHR", "SEIQHRF")) { |
202 | 2x |
dat$epi[[outName_c[1]]] <- c(0, nQuar) |
203 | 2x |
dat$epi[[outName_d[1]]] <- c(0, nHosp) |
204 |
} |
|
205 | 2x |
if (type %in% c("SEIQHRF")) { |
206 | 2x |
dat$epi[[outName_e[1]]] <- c(0, nFat) |
207 |
} |
|
208 |
} else { |
|
209 | ! |
dat$epi[[outName_a[1]]][at] <- nRecov |
210 | ! |
if (type %in% c("SEIR", "SEIQHR")) { |
211 | ! |
dat$epi[[outName_b[1]]][at] <- nProg |
212 |
} |
|
213 | ! |
if (type %in% c("SEIQHR", "SEIQHRF")) { |
214 | ! |
dat$epi[[outName_c[1]]][at] <- nQuar |
215 | ! |
dat$epi[[outName_d[1]]][at] <- nHosp |
216 |
} |
|
217 | ! |
if (type %in% c("SEIQHRF")) { |
218 | ! |
dat$epi[[outName_e[1]]][at] <- nFat |
219 |
} |
|
220 |
} |
|
221 |
|
|
222 | 2x |
return(dat) |
223 |
} |
1 |
## Churches' original function for checking consistency |
|
2 |
## internal in sirplus |
|
3 | ||
4 |
initialize.icm <- function(param, init, control, seed = NULL) { |
|
5 |
|
|
6 | 10x |
if(!is.null(seed)) set.seed(seed) |
7 |
|
|
8 |
## Master List for Data ## |
|
9 | 10x |
dat <- list() |
10 | 10x |
dat$param <- param |
11 | 10x |
dat$init <- init |
12 | 10x |
dat$control <- control |
13 |
|
|
14 |
|
|
15 |
# Set attributes |
|
16 | 10x |
dat$attr <- list() |
17 | 10x |
numeric.init <- init[which(sapply(init, class) == "numeric")] |
18 | 10x |
n <- do.call("sum", numeric.init) |
19 | 10x |
dat$attr$active <- rep(1, n) |
20 | 10x |
if (dat$param$groups == 1) { |
21 | 10x |
dat$attr$group <- rep(1, n) |
22 |
} else { |
|
23 | ! |
g2inits <- grep(".g2", names(numeric.init)) |
24 | ! |
g1inits <- setdiff(1:length(numeric.init), g2inits) |
25 | ! |
nG1 <- sum(sapply(g1inits, function(x) init[[x]])) |
26 | ! |
nG2 <- sum(sapply(g2inits, function(x) init[[x]])) |
27 | ! |
dat$attr$group <- c(rep(1, nG1), rep(2, max(0, nG2))) |
28 |
} |
|
29 |
|
|
30 |
# Initialize status and infection time |
|
31 | 10x |
dat <- init_status.icm(dat) |
32 |
|
|
33 |
|
|
34 |
# Summary out list |
|
35 | 10x |
dat <- get_prev.icm(dat, at = 1) |
36 |
|
|
37 | 10x |
return(dat) |
38 |
} |
|
39 | ||
40 | ||
41 |
init_status.icm <- function(dat) { |
|
42 |
|
|
43 |
# Variables --------------------------------------------------------------- |
|
44 | 23x |
type <- dat$control$type |
45 | 23x |
group <- dat$attr$group |
46 | 23x |
nGroups <- dat$param$groups |
47 |
|
|
48 | 23x |
nG1 <- sum(group == 1) |
49 | 23x |
nG2 <- sum(group == 2) |
50 |
|
|
51 | 23x |
e.num <- dat$init$e.num |
52 | 23x |
i.num <- dat$init$i.num |
53 | 23x |
q.num <- dat$init$q.num |
54 | 23x |
h.num <- dat$init$h.num |
55 | 23x |
r.num <- dat$init$r.num |
56 | 23x |
f.num <- dat$init$f.num |
57 | 23x |
e.num.g2 <- dat$init$e.num.g2 |
58 | 23x |
i.num.g2 <- dat$init$i.num.g2 |
59 | 23x |
q.num.g2 <- dat$init$q.num.g2 |
60 | 23x |
h.num.g2 <- dat$init$h.num.g2 |
61 | 23x |
r.num.g2 <- dat$init$r.num.g2 |
62 | 23x |
f.num.g2 <- dat$init$f.num.g2 |
63 |
|
|
64 |
# Status ------------------------------------------------------------------ |
|
65 | 23x |
status <- rep("s", nG1 + nG2) |
66 | 23x |
status[sample(which(group == 1), size = i.num)] <- "i" |
67 | 23x |
if (nGroups == 2) { |
68 | ! |
status[sample(which(group == 2), size = i.num.g2)] <- "i" |
69 |
} |
|
70 | 23x |
if (type %in% c("SIR", "SEIR", "SEIQHR", "SEIQHRF")) { |
71 | 23x |
status[sample(which(group == 1 & status == "s"), size = r.num)] <- "r" |
72 | 23x |
if (nGroups == 2) { |
73 | ! |
status[sample(which(group == 2 & status == "s"), size = r.num.g2)] <- "r" |
74 |
} |
|
75 |
} |
|
76 | 23x |
if (type %in% c("SEIR", "SEIQHR", "SEIQHRF")) { |
77 | 23x |
status[sample(which(group == 1 & status == "s"), size = e.num)] <- "e" |
78 | 23x |
if (nGroups == 2) { |
79 | ! |
status[sample(which(group == 2 & status == "s"), size = e.num.g2)] <- "e" |
80 |
} |
|
81 |
} |
|
82 | 23x |
if (type %in% c("SEIQHR", "SEIQHRF")) { |
83 | 23x |
status[sample(which(group == 1 & status == "s"), size = q.num)] <- "q" |
84 | 23x |
if (nGroups == 2) { |
85 | ! |
status[sample(which(group == 2 & status == "s"), size = q.num.g2)] <- "q" |
86 |
} |
|
87 | 23x |
status[sample(which(group == 1 & status == "s"), size = h.num)] <- "h" |
88 | 23x |
if (nGroups == 2) { |
89 | ! |
status[sample(which(group == 2 & status == "s"), size = h.num.g2)] <- "h" |
90 |
} |
|
91 |
} |
|
92 | 23x |
if (type %in% c("SEIQHRF")) { |
93 | 23x |
status[sample(which(group == 1 & status == "s"), size = f.num)] <- "f" |
94 | 23x |
if (nGroups == 2) { |
95 | ! |
status[sample(which(group == 2 & status == "s"), size = f.num.g2)] <- "f" |
96 |
} |
|
97 |
} |
|
98 |
|
|
99 | 23x |
dat$attr$status <- status |
100 |
|
|
101 |
|
|
102 |
# Exposure Time ---------------------------------------------------------- |
|
103 | 23x |
idsExp <- which(status == "e") |
104 | 23x |
expTime <- rep(NA, length(status)) |
105 |
# leave exposure time uninitialised for now, and |
|
106 |
# just set to NA at start. |
|
107 | 23x |
dat$attr$expTime <- expTime |
108 |
|
|
109 |
# Infection Time ---------------------------------------------------------- |
|
110 | 23x |
idsInf <- which(status == "i") |
111 | 23x |
infTime <- rep(NA, length(status)) |
112 | 23x |
dat$attr$infTime <- infTime # overwritten below |
113 |
|
|
114 |
# Recovery Time ---------------------------------------------------------- |
|
115 | 23x |
idsRecov <- which(status == "r") |
116 | 23x |
recovTime <- rep(NA, length(status)) |
117 | 23x |
dat$attr$recovTime <- recovTime |
118 |
|
|
119 |
# Need for Hospitalisation Time ---------------------------------------------------------- |
|
120 | 23x |
idsHosp <- which(status == "h") |
121 | 23x |
hospTime <- rep(NA, length(status)) |
122 | 23x |
dat$attr$hospTime <- hospTime |
123 |
|
|
124 |
# Quarantine Time ---------------------------------------------------------- |
|
125 | 23x |
idsQuar <- which(status == "q") |
126 | 23x |
quarTime <- rep(NA, length(status)) |
127 | 23x |
dat$attr$quarTime <- quarTime |
128 |
|
|
129 |
# Hospital-need cessation Time ---------------------------------------------------------- |
|
130 | 23x |
dischTime <- rep(NA, length(status)) |
131 | 23x |
dat$attr$dischTime <- dischTime |
132 |
|
|
133 |
# Case-fatality Time ---------------------------------------------------------- |
|
134 | 23x |
fatTime <- rep(NA, length(status)) |
135 | 23x |
dat$attr$fatTime <- fatTime |
136 |
|
|
137 |
# If vital=TRUE, infTime is a uniform draw over the duration of infection |
|
138 |
# note the initial infections may have negative infTime! |
|
139 | 23x |
if (FALSE) { |
140 |
# not sure what the following section is trying to do, but it |
|
141 |
# mucks up the gamma-distributed incumabtion periods, so set |
|
142 |
# infTime for initial infected people to t=1 instead |
|
143 | ! |
if (dat$param$vital == TRUE && dat$param$di.rate > 0) { |
144 | ! |
infTime[idsInf] <- -rgeom(n = length(idsInf), prob = dat$param$di.rate) + 2 |
145 |
} else { |
|
146 | ! |
if (dat$control$type == "SI" || dat$param$rec.rate == 0) { |
147 |
# infTime a uniform draw over the number of sim time steps |
|
148 | ! |
infTime[idsInf] <- ssample(1:(-dat$control$nsteps + 2), |
149 | ! |
length(idsInf), replace = TRUE) |
150 |
} else { |
|
151 | ! |
if (nGroups == 1) { |
152 | ! |
infTime[idsInf] <- ssample(1:(-round(1 / dat$param$rec.rate) + 2), |
153 | ! |
length(idsInf), replace = TRUE) |
154 |
} |
|
155 | ! |
if (nGroups == 2) { |
156 | ! |
infG1 <- which(status == "i" & group == 1) |
157 | ! |
infTime[infG1] <- ssample(1:(-round(1 / dat$param$rec.rate) + 2), |
158 | ! |
length(infG1), replace = TRUE) |
159 | ! |
infG2 <- which(status == "i" & group == 2) |
160 | ! |
infTime[infG2] <- ssample(1:(-round(1 / dat$param$rec.rate.g2) + 2), |
161 | ! |
length(infG2), replace = TRUE) |
162 |
} |
|
163 |
} |
|
164 |
} |
|
165 |
} |
|
166 | 23x |
infTime[idsInf] <- 1 |
167 | 23x |
dat$attr$infTime <- infTime |
168 |
|
|
169 | 23x |
return(dat) |
170 |
} |
|
171 | ||
172 | ||
173 |
infection.seiqhrf.icm <- function(dat, at, seed = NULL) { |
|
174 |
|
|
175 | 10x |
if(!is.null(seed)) set.seed(seed) |
176 |
|
|
177 | 10x |
type <- dat$control$type |
178 |
|
|
179 |
# the following checks need to be moved to control.icm in due course |
|
180 | 10x |
nsteps <- dat$control$nsteps |
181 |
|
|
182 | 10x |
act.rate.i <- dat$param$act.rate.i |
183 | 10x |
if (!(length(act.rate.i) == 1 || length(act.rate.i == nsteps))) { |
184 | ! |
stop("Length of act.rate.i must be 1 or the value of nsteps") |
185 |
} |
|
186 | 10x |
act.rate.i.g2 <- dat$param$act.rate.i.g2 |
187 | 10x |
if (!is.null(act.rate.i.g2) && |
188 | 10x |
!(length(act.rate.i.g2) == 1 || length(act.rate.i.g2 == nsteps))) { |
189 | ! |
stop("Length of act.rate.i.g2 must be 1 or the value of nsteps") |
190 |
} |
|
191 | 10x |
inf.prob.i <- dat$param$inf.prob.i |
192 | 10x |
if (!(length(inf.prob.i) == 1 || length(inf.prob.i == nsteps))) { |
193 | ! |
stop("Length of inf.prob.i must be 1 or the value of nsteps") |
194 |
} |
|
195 | 10x |
inf.prob.i.g2 <- dat$param$inf.prob.i.g2 |
196 | 10x |
if (!is.null(inf.prob.i.g2) && |
197 | 10x |
!(length(inf.prob.i.g2) == 1 || length(inf.prob.i.g2 == nsteps))) { |
198 | ! |
stop("Length of inf.prob.i.g2 must be 1 or the value of nsteps") |
199 |
} |
|
200 | 10x |
if (type %in% c("SEIQHR", "SEIQHRF")) { |
201 | 10x |
quar.rate <- dat$param$quar.rate |
202 | 10x |
if (!(length(quar.rate) == 1 || length(quar.rate == nsteps))) { |
203 | ! |
stop("Length of quar.rate must be 1 or the value of nsteps") |
204 |
} |
|
205 | 10x |
quar.rate.g2 <- dat$param$quar.rate.g2 |
206 | 10x |
if (!is.null(quar.rate.g2) && |
207 | 10x |
!(length(quar.rate.g2) == 1 || length(quar.rate.g2 == nsteps))) { |
208 | ! |
stop("Length of quar.rate.g2 must be 1 or the value of nsteps") |
209 |
} |
|
210 | 10x |
disch.rate <- dat$param$disch.rate |
211 | 10x |
if (!(length(disch.rate) == 1 || length(disch.rate == nsteps))) { |
212 | ! |
stop("Length of disch.rate must be 1 or the value of nsteps") |
213 |
} |
|
214 | 10x |
disch.rate.g2 <- dat$param$disch.rate.g2 |
215 | 10x |
if (!is.null(disch.rate.g2) && |
216 | 10x |
!(length(disch.rate.g2) == 1 || length(disch.rate.g2 == nsteps))) { |
217 | ! |
stop("Length of disch.rate.g2 must be 1 or the value of nsteps") |
218 |
} |
|
219 |
} |
|
220 |
|
|
221 | 10x |
if (type %in% c("SEIQHR", "SEIQHRF")) { |
222 | 10x |
act.rate.e <- dat$param$act.rate.e |
223 | 10x |
if (!(length(act.rate.e) == 1 || length(act.rate.e == nsteps))) { |
224 | ! |
stop("Length of act.rate.e must be 1 or the value of nsteps") |
225 |
} |
|
226 | 10x |
act.rate.e.g2 <- dat$param$act.rate.e.g2 |
227 | 10x |
if (!is.null(act.rate.e.g2) && |
228 | 10x |
!(length(act.rate.e.g2) == 1 || length(act.rate.e.g2 == nsteps))) { |
229 | ! |
stop("Length of act.rate.e.g2 must be 1 or the value of nsteps") |
230 |
} |
|
231 | 10x |
inf.prob.e <- dat$param$inf.prob.e |
232 | 10x |
if (!(length(inf.prob.e) == 1 || length(inf.prob.e == nsteps))) { |
233 | ! |
stop("Length of inf.prob.e must be 1 or the value of nsteps") |
234 |
} |
|
235 | 10x |
inf.prob.e.g2 <- dat$param$inf.prob.e.g2 |
236 | 10x |
if (!is.null(inf.prob.e.g2) && |
237 | 10x |
!(length(inf.prob.e.g2) == 1 || length(inf.prob.e.g2 == nsteps))) { |
238 | ! |
stop("Length of inf.prob.e.g2 must be 1 or the value of nsteps") |
239 |
} |
|
240 |
|
|
241 | 10x |
act.rate.q <- dat$param$act.rate.q |
242 | 10x |
if (!(length(act.rate.q) == 1 || length(act.rate.q == nsteps))) { |
243 | ! |
stop("Length of act.rate.q must be 1 or the value of nsteps") |
244 |
} |
|
245 | 10x |
act.rate.q.g2 <- dat$param$act.rate.q.g2 |
246 | 10x |
if (!is.null(act.rate.q.g2) && |
247 | 10x |
!(length(act.rate.q.g2) == 1 || length(act.rate.q.g2 == nsteps))) { |
248 | ! |
stop("Length of act.rate.q.g2 must be 1 or the value of nsteps") |
249 |
} |
|
250 | 10x |
inf.prob.q <- dat$param$inf.prob.q |
251 | 10x |
if (!(length(inf.prob.q) == 1 || length(inf.prob.q == nsteps))) { |
252 | ! |
stop("Length of inf.prob.q must be 1 or the value of nsteps") |
253 |
} |
|
254 | 10x |
inf.prob.q.g2 <- dat$param$inf.prob.q.g2 |
255 | 10x |
if (!is.null(inf.prob.q.g2) && |
256 | 10x |
!(length(inf.prob.q.g2) == 1 || length(inf.prob.q.g2 == nsteps))) { |
257 | ! |
stop("Length of inf.prob.q.g2 must be 1 or the value of nsteps") |
258 |
} |
|
259 |
} |
|
260 |
|
|
261 |
# Transmission from infected |
|
262 |
## Expected acts |
|
263 | 10x |
if (dat$param$groups == 1) { |
264 | 10x |
if (length(act.rate.i) > 1) { |
265 | ! |
acts <- round(act.rate.i[at - 1] * dat$epi$num[at - 1] / 2) |
266 |
} else { |
|
267 | 10x |
acts <- round(act.rate.i * dat$epi$num[at - 1] / 2) |
268 |
} |
|
269 |
} |
|
270 | 10x |
if (dat$param$groups == 2) { |
271 | ! |
if (dat$param$balance == "g1") { |
272 | ! |
if (length(act.rate.i) > 1) { |
273 | ! |
acts <- round(act.rate.i[at - 1] * |
274 | ! |
(dat$epi$num[at - 1] + dat$epi$num.g2[at - 1]) / 2) |
275 |
} else { |
|
276 | ! |
acts <- round(act.rate.i * |
277 | ! |
(dat$epi$num[at - 1] + dat$epi$num.g2[at - 1]) / 2) |
278 |
} |
|
279 |
} |
|
280 | ! |
if (dat$param$balance == "g2") { |
281 | ! |
if (length(act.rate.i.g2) > 1) { |
282 | ! |
acts <- round(act.rate.i.g2[at - 1] * |
283 | ! |
(dat$epi$num[at - 1] + dat$epi$num.g2[at - 1]) / 2) |
284 |
} else { |
|
285 | ! |
acts <- round(act.rate.i.g2 * |
286 | ! |
(dat$epi$num[at - 1] + dat$epi$num.g2[at - 1]) / 2) |
287 |
} |
|
288 |
} |
|
289 |
} |
|
290 |
|
|
291 |
## Edgelist |
|
292 | 10x |
if (dat$param$groups == 1) { |
293 | 10x |
p1 <- ssample(which(dat$attr$active == 1 & dat$attr$status != "f"), acts, replace = TRUE) |
294 | 10x |
p2 <- ssample(which(dat$attr$active == 1 & dat$attr$status != "f"), acts, replace = TRUE) |
295 |
} else { |
|
296 | ! |
p1 <- ssample(which(dat$attr$active == 1 & dat$attr$group == 1 & dat$attr$status != "f"), |
297 | ! |
acts, replace = TRUE) |
298 | ! |
p2 <- ssample(which(dat$attr$active == 1 & dat$attr$group == 2 & dat$attr$status != "f"), |
299 | ! |
acts, replace = TRUE) |
300 |
} |
|
301 |
|
|
302 | 10x |
del <- NULL |
303 | 10x |
if (length(p1) > 0 & length(p2) > 0) { |
304 | 10x |
del <- data.frame(p1, p2) |
305 | 10x |
if (dat$param$groups == 1) { |
306 | 10x |
while (any(del$p1 == del$p2)) { |
307 | 10x |
del$p2 <- ifelse(del$p1 == del$p2, |
308 | 10x |
ssample(which(dat$attr$active == 1 & dat$attr$status != "f"), 1), del$p2) |
309 |
} |
|
310 |
} |
|
311 |
} |
|
312 |
|
|
313 |
## Discordant edgelist (del) |
|
314 | 10x |
del$p1.stat <- dat$attr$status[del$p1] |
315 | 10x |
del$p2.stat <- dat$attr$status[del$p2] |
316 | 10x |
serodis <- (del$p1.stat == "s" & del$p2.stat == "i") | |
317 | 10x |
(del$p1.stat == "i" & del$p2.stat == "s") |
318 | 10x |
del <- del[serodis == TRUE, ] |
319 |
|
|
320 |
## Transmission on edgelist |
|
321 | 10x |
if (nrow(del) > 0) { |
322 | 10x |
if (dat$param$groups == 1) { |
323 | 10x |
if (length(inf.prob.i) > 1) { |
324 | ! |
del$tprob <- inf.prob.i[at] |
325 |
} else { |
|
326 | 10x |
del$tprob <- inf.prob.i |
327 |
} |
|
328 |
} else { |
|
329 | ! |
if (length(inf.prob.i) > 1) { |
330 | ! |
del$tprob <- ifelse(del$p1.stat == "s", inf.prob.i[at], |
331 | ! |
inf.prob.i.g2[at]) |
332 |
} else { |
|
333 | ! |
del$tprob <- ifelse(del$p1.stat == "s", inf.prob.i, |
334 | ! |
inf.prob.i.g2) |
335 |
} |
|
336 |
} |
|
337 | 10x |
if (!is.null(dat$param$inter.eff.i) && at >= dat$param$inter.start.i && |
338 | 10x |
at <= dat$param$inter.stop.i) { |
339 | ! |
del$tprob <- del$tprob * (1 - dat$param$inter.eff.i) |
340 |
} |
|
341 | 10x |
del$trans <- rbinom(nrow(del), 1, del$tprob) |
342 | 10x |
del <- del[del$trans == TRUE, ] |
343 | 10x |
if (nrow(del) > 0) { |
344 | 8x |
if (dat$param$groups == 1) { |
345 | 8x |
newIds <- unique(ifelse(del$p1.stat == "s", del$p1, del$p2)) |
346 | 8x |
nExp.i <- length(newIds) |
347 |
} |
|
348 | 8x |
if (dat$param$groups == 2) { |
349 | ! |
newIdsg1 <- unique(del$p1[del$p1.stat == "s"]) |
350 | ! |
newIdsg2 <- unique(del$p2[del$p2.stat == "s"]) |
351 | ! |
nExp.i <- length(newIdsg1) |
352 | ! |
nExpg2.i <- length(newIdsg2) |
353 | ! |
newIds <- c(newIdsg1, newIdsg2) |
354 |
} |
|
355 | 8x |
dat$attr$status[newIds] <- "e" |
356 | 8x |
dat$attr$expTime[newIds] <- at |
357 |
} else { |
|
358 | 2x |
nExp.i <- nExpg2.i <- 0 |
359 |
} |
|
360 |
} else { |
|
361 | ! |
nExp.i <- nExpg2.i <- 0 |
362 |
} |
|
363 |
|
|
364 | 10x |
if (type %in% c("SEIQHRF")) { |
365 |
|
|
366 |
# Transmission from exposed |
|
367 |
## Expected acts |
|
368 | 10x |
if (dat$param$groups == 1) { |
369 | 10x |
if (length(act.rate.e) > 1) { |
370 | ! |
acts <- round(act.rate.e[at - 1] * dat$epi$num[at - 1] / 2) |
371 |
} else { |
|
372 | 10x |
acts <- round(act.rate.e * dat$epi$num[at - 1] / 2) |
373 |
} |
|
374 |
} |
|
375 | 10x |
if (dat$param$groups == 2) { |
376 | ! |
if (dat$param$balance == "g1") { |
377 | ! |
if (length(act.rate.e.g2) > 1) { |
378 | ! |
acts <- round(act.rate.e.g2[at - 1] * |
379 | ! |
(dat$epi$num[at - 1] + dat$epi$num.g2[at - 1]) / 2) |
380 |
} else { |
|
381 | ! |
acts <- round(act.rate.e.g2 * |
382 | ! |
(dat$epi$num[at - 1] + dat$epi$num.g2[at - 1]) / 2) |
383 |
} |
|
384 |
} |
|
385 | ! |
if (dat$param$balance == "g2") { |
386 | ! |
if (length(act.rate.e.g2) > 1) { |
387 | ! |
acts <- round(act.rate.e.g2[at - 1] * |
388 | ! |
(dat$epi$num[at - 1] + dat$epi$num.g2[at - 1]) / 2) |
389 |
} else { |
|
390 | ! |
acts <- round(act.rate.e.g2 * |
391 | ! |
(dat$epi$num[at - 1] + dat$epi$num.g2[at - 1]) / 2) |
392 |
} |
|
393 |
} |
|
394 |
} |
|
395 |
|
|
396 |
## Edgelist |
|
397 | 10x |
if (dat$param$groups == 1) { |
398 | 10x |
p1 <- ssample(which(dat$attr$active == 1 & dat$attr$status != "f"), acts, replace = TRUE) |
399 | 10x |
p2 <- ssample(which(dat$attr$active == 1 & dat$attr$status != "f"), acts, replace = TRUE) |
400 |
} else { |
|
401 | ! |
p1 <- ssample(which(dat$attr$active == 1 & dat$attr$group == 1 & dat$attr$status != "f"), |
402 | ! |
acts, replace = TRUE) |
403 | ! |
p2 <- ssample(which(dat$attr$active == 1 & dat$attr$group == 2 & dat$attr$status != "f"), |
404 | ! |
acts, replace = TRUE) |
405 |
} |
|
406 |
|
|
407 | 10x |
del <- NULL |
408 | 10x |
if (length(p1) > 0 & length(p2) > 0) { |
409 | 10x |
del <- data.frame(p1, p2) |
410 | 10x |
if (dat$param$groups == 1) { |
411 | 10x |
while (any(del$p1 == del$p2)) { |
412 | 10x |
del$p2 <- ifelse(del$p1 == del$p2, |
413 | 10x |
ssample(which(dat$attr$active == 1 & dat$attr$status != "f"), 1), del$p2) |
414 |
} |
|
415 |
} |
|
416 |
|
|
417 |
## Discordant edgelist (del) |
|
418 | 10x |
del$p1.stat <- dat$attr$status[del$p1] |
419 | 10x |
del$p2.stat <- dat$attr$status[del$p2] |
420 |
# serodiscordance |
|
421 | 10x |
serodis <- (del$p1.stat == "s" & del$p2.stat == "e") | |
422 | 10x |
(del$p1.stat == "e" & del$p2.stat == "s") |
423 | 10x |
del <- del[serodis == TRUE, ] |
424 |
|
|
425 |
## Transmission on edgelist |
|
426 | 10x |
if (nrow(del) > 0) { |
427 | 8x |
if (dat$param$groups == 1) { |
428 | 8x |
if (length(inf.prob.e) > 1) { |
429 | ! |
del$tprob <- inf.prob.e[at] |
430 |
} else { |
|
431 | 8x |
del$tprob <- inf.prob.e |
432 |
} |
|
433 |
} else { |
|
434 | ! |
if (length(inf.prob.e) > 1) { |
435 | ! |
del$tprob <- ifelse(del$p1.stat == "s", inf.prob.e[at], |
436 | ! |
inf.prob.e.g2[at]) |
437 |
} else { |
|
438 | ! |
del$tprob <- ifelse(del$p1.stat == "s", inf.prob.e, |
439 | ! |
inf.prob.e.g2) |
440 |
} |
|
441 |
} |
|
442 | 8x |
if (!is.null(dat$param$inter.eff.e) && at >= dat$param$inter.start.e && |
443 | 8x |
at <= dat$param$inter.stop.e) { |
444 | ! |
del$tprob <- del$tprob * (1 - dat$param$inter.eff.e) |
445 |
} |
|
446 | 8x |
del$trans <- rbinom(nrow(del), 1, del$tprob) |
447 | 8x |
del <- del[del$trans == TRUE, ] |
448 | 8x |
if (nrow(del) > 0) { |
449 | 4x |
if (dat$param$groups == 1) { |
450 | 4x |
newIds <- unique(ifelse(del$p1.stat == "s", del$p1, del$p2)) |
451 | 4x |
nExp.e <- length(newIds) |
452 |
} |
|
453 | 4x |
if (dat$param$groups == 2) { |
454 | ! |
newIdsg1 <- unique(del$p1[del$p1.stat == "s"]) |
455 | ! |
newIdsg2 <- unique(del$p2[del$p2.stat == "s"]) |
456 | ! |
nExp.e <- length(newIdsg1) |
457 | ! |
nExpg2.e <- length(newIdsg2) |
458 | ! |
newIds <- c(newIdsg1, newIdsg2) |
459 |
} |
|
460 | 4x |
dat$attr$status[newIds] <- "e" |
461 | 4x |
dat$attr$expTime[newIds] <- at |
462 |
} else { |
|
463 | 4x |
nExp.e <- nExpg2.e <- 0 |
464 |
} |
|
465 |
} else { |
|
466 | 2x |
nExp.e <- nExpg2.e <- 0 |
467 |
} |
|
468 |
} else { |
|
469 | ! |
nExp.e <- nExpg2.e <- 0 |
470 |
} |
|
471 |
|
|
472 |
|
|
473 |
# Transmission from quarantined |
|
474 |
## Expected acts |
|
475 | 10x |
if (dat$param$groups == 1) { |
476 | 10x |
if (length(act.rate.q) > 1) { |
477 | ! |
acts <- round(act.rate.q[at - 1] * dat$epi$num[at - 1] / 2) |
478 |
} else { |
|
479 | 10x |
acts <- round(act.rate.q * dat$epi$num[at - 1] / 2) |
480 |
} |
|
481 |
} |
|
482 | 10x |
if (dat$param$groups == 2) { |
483 | ! |
if (dat$param$balance == "g1") { |
484 | ! |
if (length(act.rate.q.g2) > 1) { |
485 | ! |
acts <- round(act.rate.q.g2[at - 1] * |
486 | ! |
(dat$epi$num[at - 1] + dat$epi$num.g2[at - 1]) / 2) |
487 |
} else { |
|
488 | ! |
acts <- round(act.rate.q.g2 * |
489 | ! |
(dat$epi$num[at - 1] + dat$epi$num.g2[at - 1]) / 2) |
490 |
} |
|
491 |
} |
|
492 | ! |
if (dat$param$balance == "g2") { |
493 | ! |
if (length(act.rate.q.g2) > 1) { |
494 | ! |
acts <- round(act.rate.q.g2[at - 1] * |
495 | ! |
(dat$epi$num[at - 1] + dat$epi$num.g2[at - 1]) / 2) |
496 |
} else { |
|
497 | ! |
acts <- round(act.rate.q.g2 * |
498 | ! |
(dat$epi$num[at - 1] + dat$epi$num.g2[at - 1]) / 2) |
499 |
} |
|
500 |
} |
|
501 |
} |
|
502 |
|
|
503 |
## Edgelist |
|
504 | 10x |
if (dat$param$groups == 1) { |
505 | 10x |
p1 <- ssample(which(dat$attr$active == 1 & dat$attr$status != "f"), acts, replace = TRUE) |
506 | 10x |
p2 <- ssample(which(dat$attr$active == 1 & dat$attr$status != "f"), acts, replace = TRUE) |
507 |
} else { |
|
508 | ! |
p1 <- ssample(which(dat$attr$active == 1 & dat$attr$group == 1 & dat$attr$status != "f"), |
509 | ! |
acts, replace = TRUE) |
510 | ! |
p2 <- ssample(which(dat$attr$active == 1 & dat$attr$group == 2 & dat$attr$status != "f"), |
511 | ! |
acts, replace = TRUE) |
512 |
} |
|
513 |
|
|
514 | 10x |
del <- NULL |
515 | 10x |
if (length(p1) > 0 & length(p2) > 0) { |
516 | 10x |
del <- data.frame(p1, p2) |
517 | 10x |
if (dat$param$groups == 1) { |
518 | 10x |
while (any(del$p1 == del$p2)) { |
519 | 6x |
del$p2 <- ifelse(del$p1 == del$p2, |
520 | 6x |
ssample(which(dat$attr$active == 1 & dat$attr$status != "f"), 1), del$p2) |
521 |
} |
|
522 |
} |
|
523 |
|
|
524 |
## Discordant edgelist (del) |
|
525 | 10x |
del$p1.stat <- dat$attr$status[del$p1] |
526 | 10x |
del$p2.stat <- dat$attr$status[del$p2] |
527 |
# serodiscordance |
|
528 | 10x |
serodis <- (del$p1.stat == "s" & del$p2.stat == "q") | |
529 | 10x |
(del$p1.stat == "q" & del$p2.stat == "s") |
530 | 10x |
del <- del[serodis == TRUE, ] |
531 |
|
|
532 |
## Transmission on edgelist |
|
533 | 10x |
if (nrow(del) > 0) { |
534 | ! |
if (dat$param$groups == 1) { |
535 | ! |
if (length(inf.prob.q) > 1) { |
536 | ! |
del$tprob <- inf.prob.q[at] |
537 |
} else { |
|
538 | ! |
del$tprob <- inf.prob.q |
539 |
} |
|
540 |
} else { |
|
541 | ! |
if (length(inf.prob.q) > 1) { |
542 | ! |
del$tprob <- ifelse(del$p1.stat == "s", inf.prob.q[at], |
543 | ! |
inf.prob.q.g2[at]) |
544 |
} else { |
|
545 | ! |
del$tprob <- ifelse(del$p1.stat == "s", inf.prob.q, |
546 | ! |
inf.prob.q.g2) |
547 |
} |
|
548 |
} |
|
549 | ! |
if (!is.null(dat$param$inter.eff.q) && at >= dat$param$inter.start.q && |
550 | ! |
at <= dat$param$inter.stop.q) { |
551 | ! |
del$tprob <- del$tprob * (1 - dat$param$inter.eff.q) |
552 |
} |
|
553 | ! |
del$trans <- rbinom(nrow(del), 1, del$tprob) |
554 | ! |
del <- del[del$trans == TRUE, ] |
555 | ! |
if (nrow(del) > 0) { |
556 | ! |
if (dat$param$groups == 1) { |
557 | ! |
newIds <- unique(ifelse(del$p1.stat == "s", del$p1, del$p2)) |
558 | ! |
nExp.q <- length(newIds) |
559 |
} |
|
560 | ! |
if (dat$param$groups == 2) { |
561 | ! |
newIdsg1 <- unique(del$p1[del$p1.stat == "s"]) |
562 | ! |
newIdsg2 <- unique(del$p2[del$p2.stat == "s"]) |
563 | ! |
nExp.q <- length(newIdsg1) |
564 | ! |
nExpg2.q <- length(newIdsg2) |
565 | ! |
newIds <- c(newIdsg1, newIdsg2) |
566 |
} |
|
567 | ! |
dat$attr$status[newIds] <- "e" |
568 | ! |
dat$attr$expTime[newIds] <- at |
569 |
} else { |
|
570 | ! |
nExp.q <- nExpg2.q <- 0 |
571 |
} |
|
572 |
} else { |
|
573 | 10x |
nExp.q <- nExpg2.q <- 0 |
574 |
} |
|
575 |
} else { |
|
576 | ! |
nExp.q <- nExpg2.q <- 0 |
577 |
} |
|
578 |
} |
|
579 |
|
|
580 |
## Output |
|
581 | 10x |
if (type %in% c("SEIQHR", "SEIQHRF")) { |
582 | 10x |
if (at == 2) { |
583 | 10x |
dat$epi$se.flow <- c(0, nExp.i + nExp.q) |
584 |
} else { |
|
585 | ! |
dat$epi$se.flow[at] <- nExp.i + nExp.q |
586 |
} |
|
587 | 10x |
if (dat$param$groups == 2) { |
588 | ! |
if (at == 2) { |
589 | ! |
dat$epi$se.flow.g2 <- c(0, nExpg2.i + nExpg2.q ) |
590 |
} else { |
|
591 | ! |
dat$epi$se.flow.g2[at] <- nExpg2.i + nExpg2.q |
592 |
} |
|
593 |
} |
|
594 |
} else { |
|
595 | ! |
if (at == 2) { |
596 | ! |
dat$epi$se.flow <- c(0, nExp.i) |
597 |
} else { |
|
598 | ! |
dat$epi$se.flow[at] <- nExp.i |
599 |
} |
|
600 | ! |
if (dat$param$groups == 2) { |
601 | ! |
if (at == 2) { |
602 | ! |
dat$epi$se.flow.g2 <- c(0, nExpg2.i) |
603 |
} else { |
|
604 | ! |
dat$epi$se.flow.g2[at] <- nExpg2.i |
605 |
} |
|
606 |
} |
|
607 |
} |
|
608 | 10x |
return(dat) |
609 |
|
|
610 |
} |
|
611 | ||
612 | ||
613 | ||
614 | ||
615 |
progress.seiqhrf.icm <- function(dat, at, seed = NULL) { |
|
616 |
|
|
617 | ! |
if(!is.null(seed)) set.seed(seed) |
618 |
#print(at) |
|
619 |
#print(dat$control$type) |
|
620 |
#print("-------") |
|
621 |
|
|
622 |
# Conditions -------------------------------------------------------------- |
|
623 | ! |
if (!(dat$control$type %in% c("SIR", "SIS", "SEIR", "SEIQHR", "SEIQHRF"))) { |
624 | ! |
return(dat) |
625 |
} |
|
626 |
|
|
627 |
|
|
628 |
# Variables --------------------------------------------------------------- |
|
629 | ! |
active <- dat$attr$active |
630 | ! |
status <- dat$attr$status |
631 |
|
|
632 | ! |
groups <- dat$param$groups |
633 | ! |
group <- dat$attr$group |
634 |
|
|
635 | ! |
type <- dat$control$type |
636 | ! |
recovState <- ifelse(type %in% c("SIR", "SEIR", "SEIQHR", "SEIQHRF"), "r", "s") |
637 | ! |
progState <- "i" |
638 | ! |
quarState <- "q" |
639 | ! |
hospState <- "h" |
640 | ! |
fatState <- "f" |
641 |
|
|
642 |
# --- progress from exposed to infectious ---- |
|
643 | ! |
prog.rand <- dat$control$prog.rand |
644 | ! |
prog.rate <- dat$param$prog.rate |
645 | ! |
prog.rate.g2 <- dat$param$prog.rate.g2 |
646 | ! |
prog.dist.scale <- dat$param$prog.dist.scale |
647 | ! |
prog.dist.shape <- dat$param$prog.dist.shape |
648 | ! |
prog.dist.scale.g2 <- dat$param$prog.dist.scale.g2 |
649 | ! |
prog.dist.shape.g2 <- dat$param$prog.dist.shape.g2 |
650 |
|
|
651 | ! |
nProg <- nProgG2 <- 0 |
652 | ! |
idsElig <- which(active == 1 & status == "e") |
653 | ! |
nElig <- length(idsElig) |
654 |
|
|
655 | ! |
if (nElig > 0) { |
656 |
|
|
657 | ! |
gElig <- group[idsElig] |
658 | ! |
rates <- c(prog.rate, prog.rate.g2) |
659 | ! |
ratesElig <- rates[gElig] |
660 |
|
|
661 | ! |
if (prog.rand == TRUE) { |
662 | ! |
vecProg <- which(rbinom(nElig, 1, ratesElig) == 1) |
663 | ! |
if (length(vecProg) > 0) { |
664 | ! |
idsProg <- idsElig[vecProg] |
665 | ! |
nProg <- sum(group[idsProg] == 1) |
666 | ! |
nProgG2 <- sum(group[idsProg] == 2) |
667 | ! |
status[idsProg] <- progState |
668 | ! |
dat$attr$infTime[idsProg] <- at |
669 |
} |
|
670 |
} else { |
|
671 | ! |
vecTimeSinceExp <- at - dat$attr$expTime[idsElig] |
672 | ! |
gammaRatesElig <- pweibull(vecTimeSinceExp, prog.dist.shape, scale=prog.dist.scale) |
673 | ! |
nProg <- round(sum(gammaRatesElig[gElig == 1], na.rm=TRUE)) |
674 | ! |
if (nProg > 0) { |
675 | ! |
ids2bProg <- ssample(idsElig[gElig == 1], |
676 | ! |
nProg, prob = gammaRatesElig[gElig == 1]) |
677 | ! |
status[ids2bProg] <- progState |
678 | ! |
dat$attr$infTime[ids2bProg] <- at |
679 |
# debug |
|
680 | ! |
if (FALSE & at <= 30) { |
681 | ! |
print(paste("at:", at)) |
682 | ! |
print("idsElig:") |
683 | ! |
print(idsElig[gElig == 1]) |
684 | ! |
print("vecTimeSinceExp:") |
685 | ! |
print(vecTimeSinceExp[gElig == 1]) |
686 | ! |
print("gammaRatesElig:") |
687 | ! |
print(gammaRatesElig) |
688 | ! |
print(paste("nProg:",nProg)) |
689 | ! |
print(paste("sum of elig rates:", round(sum(gammaRatesElig[gElig == 1])))) |
690 | ! |
print(paste("sum(gElig == 1):", sum(gElig == 1))) |
691 | ! |
print("ids progressed:") |
692 | ! |
print(ids2bProg) |
693 | ! |
print("probs of ids to be progressed:") |
694 | ! |
print(gammaRatesElig[which(idsElig %in% ids2bProg)]) |
695 | ! |
print("days since exposed of ids to be progressed:") |
696 | ! |
print(vecTimeSinceExp[which(idsElig %in% ids2bProg)]) |
697 | ! |
print("------") |
698 |
} |
|
699 |
} |
|
700 | ! |
if (groups == 2) { |
701 | ! |
nProgG2 <- round(sum(gammaRatesElig[gElig == 2], na.rm=TRUE)) |
702 | ! |
if (nProgG2 > 0) { |
703 | ! |
ids2bProgG2 <- ssample(idsElig[gElig == 2], |
704 | ! |
nProgG2, prob = gammaRatesElig[gElig == 2]) |
705 | ! |
status[ids2bProgG2] <- progState |
706 | ! |
dat$attr$infTime[ids2bProgG2] <- at |
707 |
} |
|
708 |
} |
|
709 |
} |
|
710 |
} |
|
711 | ! |
dat$attr$status <- status |
712 |
|
|
713 | ! |
if (type %in% c("SEIQHR", "SEIQHRF")) { |
714 |
# ----- quarantine ------- |
|
715 | ! |
quar.rand <- dat$control$quar.rand |
716 | ! |
quar.rate <- dat$param$quar.rate |
717 | ! |
quar.rate.g2 <- dat$param$quar.rate.g2 |
718 |
|
|
719 | ! |
nQuar <- nQuarG2 <- 0 |
720 | ! |
idsElig <- which(active == 1 & status == "i") |
721 | ! |
nElig <- length(idsElig) |
722 |
|
|
723 | ! |
if (nElig > 0) { |
724 |
|
|
725 | ! |
gElig <- group[idsElig] |
726 | ! |
rates <- c(quar.rate, quar.rate.g2) |
727 |
|
|
728 | ! |
if (length(quar.rate) > 1) { |
729 | ! |
qrate <- quar.rate[at] |
730 |
} else { |
|
731 | ! |
qrate <- quar.rate |
732 |
} |
|
733 | ! |
if (length(quar.rate.g2) > 1) { |
734 | ! |
qrate.g2 <- quar.rate.g2[at] |
735 |
} else { |
|
736 | ! |
qrate.g2 <- quar.rate.g2 |
737 |
} |
|
738 | ! |
rates <- c(qrate, qrate.g2) |
739 | ! |
ratesElig <- rates[gElig] |
740 | ! |
if (quar.rand == TRUE) { |
741 | ! |
vecQuar <- which(rbinom(nElig, 1, ratesElig) == 1) |
742 | ! |
if (length(vecQuar) > 0) { |
743 | ! |
idsQuar <- idsElig[vecQuar] |
744 | ! |
nQuar <- sum(group[idsQuar] == 1) |
745 | ! |
nQuarG2 <- sum(group[idsQuar] == 2) |
746 | ! |
status[idsQuar] <- quarState |
747 | ! |
dat$attr$quarTime[idsQuar] <- at |
748 |
} |
|
749 |
} else { |
|
750 | ! |
nQuar <- min(round(sum(ratesElig[gElig == 1])), sum(gElig == 1)) |
751 | ! |
idsQuar <- ssample(idsElig[gElig == 1], nQuar) |
752 | ! |
status[idsQuar] <- quarState |
753 | ! |
dat$attr$quarTime[idsQuar] <- at |
754 | ! |
if (groups == 2) { |
755 | ! |
nQuarG2 <- min(round(sum(ratesElig[gElig == 2])), sum(gElig == 2)) |
756 | ! |
idsQuarG2 <- ssample(idsElig[gElig == 2], nQuarG2) |
757 | ! |
status[idsQuarG2] <- quarState |
758 | ! |
dat$attr$quarTime[idsQuarG2] <- at |
759 |
} |
|
760 |
} |
|
761 |
} |
|
762 | ! |
dat$attr$status <- status |
763 |
|
|
764 |
# ----- need to be hospitalised ------- |
|
765 | ! |
hosp.rand <- dat$control$hosp.rand |
766 | ! |
hosp.rate <- dat$param$hosp.rate |
767 | ! |
hosp.rate.g2 <- dat$param$hosp.rate.g2 |
768 |
|
|
769 | ! |
nHosp <- nHospG2 <- 0 |
770 | ! |
idsElig <- which(active == 1 & (status == "i" | status == "q")) |
771 | ! |
nElig <- length(idsElig) |
772 | ! |
idsHosp <- numeric(0) |
773 |
|
|
774 | ! |
if (nElig > 0) { |
775 |
|
|
776 | ! |
gElig <- group[idsElig] |
777 | ! |
rates <- c(hosp.rate, hosp.rate.g2) |
778 | ! |
ratesElig <- rates[gElig] |
779 |
|
|
780 | ! |
if (hosp.rand == TRUE) { |
781 | ! |
vecHosp <- which(rbinom(nElig, 1, ratesElig) == 1) |
782 | ! |
if (length(vecHosp) > 0) { |
783 | ! |
idsHosp <- idsElig[vecHosp] |
784 | ! |
nHosp <- sum(group[idsHosp] == 1) |
785 | ! |
nHospG2 <- sum(group[idsHosp] == 2) |
786 | ! |
status[idsHosp] <- hospState |
787 |
} |
|
788 |
} else { |
|
789 | ! |
nHosp <- min(round(sum(ratesElig[gElig == 1])), sum(gElig == 1)) |
790 | ! |
idsHosp <- ssample(idsElig[gElig == 1], nHosp) |
791 | ! |
status[idsHosp] <- hospState |
792 | ! |
if (groups == 2) { |
793 | ! |
nHospG2 <- min(round(sum(ratesElig[gElig == 2])), sum(gElig == 2)) |
794 | ! |
idsHospG2 <- ssample(idsElig[gElig == 2], nHospG2) |
795 | ! |
status[idsHospG2] <- hospState |
796 | ! |
idsHosp <- c(idsHosp, idsHospG2) |
797 |
} |
|
798 |
} |
|
799 |
} |
|
800 | ! |
dat$attr$status <- status |
801 | ! |
dat$attr$hospTime[idsHosp] <- at |
802 |
|
|
803 |
# ----- discharge from need to be hospitalised ------- |
|
804 | ! |
disch.rand <- dat$control$disch.rand |
805 | ! |
disch.rate <- dat$param$disch.rate |
806 | ! |
disch.rate.g2 <- dat$param$disch.rate.g2 |
807 |
|
|
808 | ! |
nDisch <- nDischG2 <- 0 |
809 | ! |
idsElig <- which(active == 1 & status == "h") |
810 | ! |
nElig <- length(idsElig) |
811 | ! |
idsDisch <- numeric(0) |
812 |
|
|
813 | ! |
if (nElig > 0) { |
814 |
|
|
815 | ! |
gElig <- group[idsElig] |
816 | ! |
rates <- c(disch.rate, disch.rate.g2) |
817 |
|
|
818 | ! |
if (length(disch.rate) > 1) { |
819 | ! |
dcrate <- disch.rate[at] |
820 |
} else { |
|
821 | ! |
dcrate <- disch.rate |
822 |
} |
|
823 | ! |
if (length(disch.rate.g2) > 1) { |
824 | ! |
dcrate.g2 <- disch.rate.g2[at] |
825 |
} else { |
|
826 | ! |
dcrate.g2 <- disch.rate.g2 |
827 |
} |
|
828 |
|
|
829 | ! |
rates <- c(dcrate, dcrate.g2) |
830 | ! |
ratesElig <- rates[gElig] |
831 |
|
|
832 | ! |
if (disch.rand == TRUE) { |
833 | ! |
vecDisch <- which(rbinom(nElig, 1, ratesElig) == 1) |
834 | ! |
if (length(vecDisch) > 0) { |
835 | ! |
idsDisch <- idsElig[vecDisch] |
836 | ! |
nDisch <- sum(group[idsDisch] == 1) |
837 | ! |
nDischG2 <- sum(group[idsDisch] == 2) |
838 | ! |
status[idsDisch] <- recovState |
839 |
} |
|
840 |
} else { |
|
841 | ! |
nDisch <- min(round(sum(ratesElig[gElig == 1])), sum(gElig == 1)) |
842 | ! |
idsDisch <- ssample(idsElig[gElig == 1], nDisch) |
843 | ! |
status[idsDisch] <- recovState |
844 | ! |
if (groups == 2) { |
845 | ! |
nDischG2 <- min(round(sum(ratesElig[gElig == 2])), sum(gElig == 2)) |
846 | ! |
idsDischG2 <- ssample(idsElig[gElig == 2], nDischG2) |
847 | ! |
status[idsDischG2] <- recovState |
848 | ! |
idsDisch <- c(idsDisch, idsDischG2) |
849 |
} |
|
850 |
} |
|
851 |
} |
|
852 | ! |
dat$attr$status <- status |
853 | ! |
dat$attr$dischTime[idsDisch] <- at |
854 |
} |
|
855 |
|
|
856 |
# ----- recover ------- |
|
857 | ! |
rec.rand <- dat$control$rec.rand |
858 | ! |
rec.rate <- dat$param$rec.rate |
859 | ! |
rec.rate.g2 <- dat$param$rec.rate.g2 |
860 | ! |
rec.dist.scale <- dat$param$rec.dist.scale |
861 | ! |
rec.dist.shape <- dat$param$rec.dist.shape |
862 | ! |
rec.dist.scale.g2 <- dat$param$rec.dist.scale.g2 |
863 | ! |
rec.dist.shape.g2 <- dat$param$rec.dist.shape.g2 |
864 |
|
|
865 | ! |
nRecov <- nRecovG2 <- 0 |
866 | ! |
idsElig <- which(active == 1 & (status == "i" | status == "q" | status == "h")) |
867 | ! |
nElig <- length(idsElig) |
868 | ! |
idsRecov <- numeric(0) |
869 |
|
|
870 | ! |
if (nElig > 0) { |
871 |
|
|
872 | ! |
gElig <- group[idsElig] |
873 | ! |
rates <- c(rec.rate, rec.rate.g2) |
874 | ! |
ratesElig <- rates[gElig] |
875 |
|
|
876 | ! |
if (rec.rand == TRUE) { |
877 | ! |
vecRecov <- which(rbinom(nElig, 1, ratesElig) == 1) |
878 | ! |
if (length(vecRecov) > 0) { |
879 | ! |
idsRecov <- idsElig[vecRecov] |
880 | ! |
nRecov <- sum(group[idsRecov] == 1) |
881 | ! |
nRecovG2 <- sum(group[idsRecov] == 2) |
882 | ! |
status[idsRecov] <- recovState |
883 |
} |
|
884 |
} else { |
|
885 | ! |
vecTimeSinceExp <- at - dat$attr$expTime[idsElig] |
886 | ! |
vecTimeSinceExp[is.na(vecTimeSinceExp)] <- 0 |
887 | ! |
gammaRatesElig <- pweibull(vecTimeSinceExp, rec.dist.shape, scale=rec.dist.scale) |
888 | ! |
nRecov <- round(sum(gammaRatesElig[gElig == 1], na.rm=TRUE)) |
889 | ! |
if (nRecov > 0) { |
890 | ! |
idsRecov <- ssample(idsElig[gElig == 1], |
891 | ! |
nRecov, prob = gammaRatesElig[gElig == 1]) |
892 | ! |
status[idsRecov] <- recovState |
893 |
# debug |
|
894 | ! |
if (FALSE & at <= 30) { |
895 | ! |
print(paste("at:", at)) |
896 | ! |
print("idsElig:") |
897 | ! |
print(idsElig[gElig == 1]) |
898 | ! |
print("vecTimeSinceExp:") |
899 | ! |
print(vecTimeSinceExp[gElig == 1]) |
900 | ! |
print("gammaRatesElig:") |
901 | ! |
print(gammaRatesElig) |
902 | ! |
print(paste("nRecov:",nRecov)) |
903 | ! |
print(paste("sum of elig rates:", round(sum(gammaRatesElig[gElig == 1])))) |
904 | ! |
print(paste("sum(gElig == 1):", sum(gElig == 1))) |
905 | ! |
print("ids recovered:") |
906 | ! |
print(idsRecov) |
907 | ! |
print("probs of ids to be progressed:") |
908 | ! |
print(gammaRatesElig[which(idsElig %in% idsRecov)]) |
909 | ! |
print("days since exposed of ids to be Recovered:") |
910 | ! |
print(vecTimeSinceExp[which(idsElig %in% idsRecov)]) |
911 | ! |
print("------") |
912 |
} |
|
913 |
|
|
914 |
} |
|
915 | ! |
if (groups == 2) { |
916 | ! |
nRecovG2 <- round(sum(gammaRatesElig[gElig == 2], na.rm=TRUE)) |
917 | ! |
if (nRecovG2 > 0) { |
918 | ! |
idsRecovG2 <- ssample(idsElig[gElig == 2], |
919 | ! |
nRecovG2, prob = gammaRatesElig[gElig == 2]) |
920 | ! |
status[idsRecovG2] <- recovState |
921 | ! |
idsRecov <- c(idsRecov, idsRecovG2) |
922 |
} |
|
923 |
} |
|
924 |
} |
|
925 |
} |
|
926 | ! |
dat$attr$status <- status |
927 | ! |
dat$attr$recovTime[idsRecov] <- at |
928 |
|
|
929 | ! |
fatEnable <- TRUE |
930 | ! |
if (fatEnable & type %in% c("SEIQHRF")) { |
931 |
# ----- case fatality ------- |
|
932 | ! |
fat.rand <- dat$control$fat.rand |
933 | ! |
fat.rate.base <- dat$param$fat.rate.base |
934 | ! |
fat.rate.base.g2 <- dat$param$fat.rate.base.g2 |
935 | ! |
fat.rate.base.g2 <- ifelse(is.null(fat.rate.base.g2), |
936 | ! |
0, fat.rate.base.g2) |
937 | ! |
fat.rate.overcap <- dat$param$fat.rate.overcap |
938 | ! |
fat.rate.overcap.g2 <- dat$param$fat.rate.overcap.g2 |
939 | ! |
fat.rate.overcap.g2 <- ifelse(is.null(fat.rate.overcap.g2), |
940 | ! |
0, fat.rate.overcap.g2) |
941 | ! |
hosp.cap <- dat$param$hosp.cap |
942 | ! |
fat.tcoeff <- dat$param$fat.tcoeff |
943 |
|
|
944 | ! |
nFat <- nFatG2 <- 0 |
945 | ! |
idsElig <- which(active == 1 & status =="h") |
946 | ! |
nElig <- length(idsElig) |
947 |
|
|
948 | ! |
if (nElig > 0) { |
949 | ! |
gElig <- group[idsElig] |
950 | ! |
timeInHospElig <- at - dat$attr$hospTime[idsElig] |
951 | ! |
rates <- c(fat.rate.base, fat.rate.base.g2) |
952 | ! |
h.num.yesterday <- 0 |
953 | ! |
if (!is.null(dat$epi$h.num[at - 1])) { |
954 | ! |
h.num.yesterday <- dat$epi$h.num[at - 1] |
955 | ! |
if (h.num.yesterday > hosp.cap) { |
956 | ! |
blended.rate <- ((hosp.cap * fat.rate.base) + |
957 | ! |
((h.num.yesterday - hosp.cap) * fat.rate.overcap)) / |
958 | ! |
h.num.yesterday |
959 | ! |
blended.rate.g2 <- ((hosp.cap * fat.rate.base.g2) + |
960 | ! |
((h.num.yesterday - hosp.cap) * fat.rate.overcap.g2)) / |
961 | ! |
h.num.yesterday |
962 | ! |
rates <- c(blended.rate, blended.rate.g2) |
963 |
} |
|
964 |
} |
|
965 | ! |
ratesElig <- rates[gElig] |
966 | ! |
ratesElig <- ratesElig + timeInHospElig*fat.tcoeff*ratesElig |
967 |
|
|
968 | ! |
if (fat.rand == TRUE) { |
969 | ! |
vecFat <- which(rbinom(nElig, 1, ratesElig) == 1) |
970 | ! |
if (length(vecFat) > 0) { |
971 | ! |
idsFat <- idsElig[vecFat] |
972 | ! |
nFat <- sum(group[idsFat] == 1) |
973 | ! |
nFatG2 <- sum(group[idsFat] == 2) |
974 | ! |
status[idsFat] <- fatState |
975 | ! |
dat$attr$fatTime[idsFat] <- at |
976 |
} |
|
977 |
} else { |
|
978 | ! |
nFat <- min(round(sum(ratesElig[gElig == 1])), sum(gElig == 1)) |
979 | ! |
idsFat <- ssample(idsElig[gElig == 1], nFat) |
980 | ! |
status[idsFat] <- fatState |
981 | ! |
dat$attr$fatTime[idsFat] <- at |
982 | ! |
if (groups == 2) { |
983 | ! |
nFatG2 <- min(round(sum(ratesElig[gElig == 2])), sum(gElig == 2)) |
984 | ! |
idsFatG2 <- ssample(idsElig[gElig == 2], nFatG2) |
985 | ! |
status[idsFatG2] <- fatState |
986 | ! |
dat$attr$fatTime[idsFatG2] <- at |
987 |
} |
|
988 |
} |
|
989 |
} |
|
990 | ! |
dat$attr$status <- status |
991 |
} |
|
992 |
|
|
993 |
# Output ------------------------------------------------------------------ |
|
994 | ! |
outName_a <- ifelse(type %in% c("SIR", "SEIR"), "ir.flow", "is.flow") |
995 | ! |
outName_a[2] <- paste0(outName_a, ".g2") |
996 | ! |
if (type %in% c("SEIR", "SEIQHR", "SEIQHRF")) { |
997 | ! |
outName_b <- "ei.flow" |
998 | ! |
outName_b[2] <- paste0(outName_b, ".g2") |
999 |
} |
|
1000 | ! |
if (type %in% c("SEIQHR", "SEIQHRF")) { |
1001 | ! |
outName_c <- "iq.flow" |
1002 | ! |
outName_c[2] <- paste0(outName_c, ".g2") |
1003 | ! |
outName_d <- "iq2h.flow" |
1004 | ! |
outName_d[2] <- paste0(outName_d, ".g2") |
1005 |
} |
|
1006 | ! |
if (type %in% c("SEIQHRF")) { |
1007 | ! |
outName_e <- "hf.flow" |
1008 | ! |
outName_e[2] <- paste0(outName_e, ".g2") |
1009 |
} |
|
1010 |
## Summary statistics |
|
1011 | ! |
if (at == 2) { |
1012 | ! |
dat$epi[[outName_a[1]]] <- c(0, nRecov) |
1013 | ! |
if (type %in% c("SEIR", "SEIQHR")) { |
1014 | ! |
dat$epi[[outName_b[1]]] <- c(0, nProg) |
1015 |
} |
|
1016 | ! |
if (type %in% c("SEIQHR", "SEIQHRF")) { |
1017 | ! |
dat$epi[[outName_c[1]]] <- c(0, nQuar) |
1018 | ! |
dat$epi[[outName_d[1]]] <- c(0, nHosp) |
1019 |
} |
|
1020 | ! |
if (fatEnable & type %in% c("SEIQHRF")) { |
1021 | ! |
dat$epi[[outName_e[1]]] <- c(0, nFat) |
1022 |
} |
|
1023 |
} else { |
|
1024 | ! |
dat$epi[[outName_a[1]]][at] <- nRecov |
1025 | ! |
if (type %in% c("SEIR", "SEIQHR")) { |
1026 | ! |
dat$epi[[outName_b[1]]][at] <- nProg |
1027 |
} |
|
1028 | ! |
if (type %in% c("SEIQHR", "SEIQHRF")) { |
1029 | ! |
dat$epi[[outName_c[1]]][at] <- nQuar |
1030 | ! |
dat$epi[[outName_d[1]]][at] <- nHosp |
1031 |
} |
|
1032 | ! |
if (fatEnable & type %in% c("SEIQHRF")) { |
1033 | ! |
dat$epi[[outName_e[1]]][at] <- nFat |
1034 |
} |
|
1035 |
} |
|
1036 | ! |
if (groups == 2) { |
1037 | ! |
if (at == 2) { |
1038 | ! |
dat$epi[[outName_a[2]]] <- c(0, nRecovG2) |
1039 | ! |
if (type %in% c("SEIR", "SEIQHR", "SEIQHRF")) { |
1040 | ! |
dat$epi[[outName_b[2]]] <- c(0, nProgG2) |
1041 |
} |
|
1042 | ! |
if (type %in% c("SEIQHR", "SEIQHRF")) { |
1043 | ! |
dat$epi[[outName_c[2]]] <- c(0, nQuarG2) |
1044 | ! |
dat$epi[[outName_d[2]]] <- c(0, nHospG2) |
1045 |
} |
|
1046 | ! |
if (type %in% c("SEIQHRF")) { |
1047 | ! |
dat$epi[[outName_e[2]]] <- c(0, nFatG2) |
1048 |
} |
|
1049 |
} else { |
|
1050 | ! |
dat$epi[[outName_a[2]]][at] <- nRecovG2 |
1051 | ! |
if (type %in% c("SEIR", "SEIQHR", "SEIQHRF")) { |
1052 | ! |
dat$epi[[outName_b[2]]][at] <- nProgG2 |
1053 |
} |
|
1054 | ! |
if (type %in% c("SEIQHR", "SEIQHRF")) { |
1055 | ! |
dat$epi[[outName_c[2]]][at] <- nQuarG2 |
1056 | ! |
dat$epi[[outName_d[2]]][at] <- nHospG2 |
1057 |
} |
|
1058 | ! |
if (type %in% c("SEIQHRF")) { |
1059 | ! |
dat$epi[[outName_e[2]]][at] <- nFatG2 |
1060 |
} |
|
1061 |
} |
|
1062 |
} |
|
1063 |
|
|
1064 | ! |
return(dat) |
1065 |
} |
|
1066 | ||
1067 | ||
1068 | ||
1069 |
arrivals.seiqhrf.icm <- function(dat, at, seed = NULL) { |
|
1070 |
|
|
1071 | 10x |
if(!is.null(seed)) set.seed(seed) |
1072 |
|
|
1073 |
# Conditions -------------------------------------------------------------- |
|
1074 | 10x |
if (dat$param$vital == FALSE) { |
1075 | ! |
return(dat) |
1076 |
} |
|
1077 |
|
|
1078 |
# Variables --------------------------------------------------------------- |
|
1079 | 10x |
a.rate <- dat$param$a.rate |
1080 | 10x |
a.rate.g2 <- dat$param$a.rate.g2 |
1081 | 10x |
a.rand <- dat$control$a.rand |
1082 | 10x |
groups <- dat$param$groups |
1083 | 10x |
nOld <- dat$epi$num[at - 1] |
1084 |
|
|
1085 |
# checking params, should be in control.icm or params.icm eventually |
|
1086 | 10x |
type <- dat$control$type |
1087 | 10x |
nsteps <- dat$control$nsteps |
1088 |
|
|
1089 | 10x |
if (!(length(a.rate) == 1 || length(a.rate == nsteps))) { |
1090 | ! |
stop("Length of a.rate must be 1 or the value of nsteps") |
1091 |
} |
|
1092 | 10x |
if (!is.null(a.rate.g2) && |
1093 | 10x |
!(length(a.rate.g2) == 1 || length(a.rate.g2 == nsteps))) { |
1094 | ! |
stop("Length of a.rate.g2 must be 1 or the value of nsteps") |
1095 |
} |
|
1096 |
|
|
1097 | 10x |
a.prop.e <- dat$param$a.prop.e |
1098 | 10x |
if (!(length(a.prop.e) == 1 || length(a.prop.e == nsteps))) { |
1099 | ! |
stop("Length of a.prop.e must be 1 or the value of nsteps") |
1100 |
} |
|
1101 | 10x |
a.prop.i <- dat$param$a.prop.i |
1102 | 10x |
if (!(length(a.prop.i) == 1 || length(a.prop.i == nsteps))) { |
1103 | ! |
stop("Length of a.prop.i must be 1 or the value of nsteps") |
1104 |
} |
|
1105 | 10x |
a.prop.q <- dat$param$a.prop.q |
1106 | 10x |
if (!(length(a.prop.q) == 1 || length(a.prop.q == nsteps))) { |
1107 | ! |
stop("Length of a.prop.q must be 1 or the value of nsteps") |
1108 |
} |
|
1109 |
|
|
1110 | 10x |
a.prop.e.g2 <- dat$param$a.prop.e.g2 |
1111 | 10x |
if (!is.null(a.prop.e.g2) && |
1112 | 10x |
!(length(a.prop.e.g2) == 1 || length(a.prop.e.g2 == nsteps))) { |
1113 | ! |
stop("Length of a.prop.e.g2 must be 1 or the value of nsteps") |
1114 |
} |
|
1115 | 10x |
a.prop.i.g2 <- dat$param$a.prop.i.g2 |
1116 | 10x |
if (!is.null(a.prop.i.g2) && |
1117 | 10x |
!(length(a.prop.i.g2) == 1 || length(a.prop.i.g2 == nsteps))) { |
1118 | ! |
stop("Length of a.prop.i.g2 must be 1 or the value of nsteps") |
1119 |
} |
|
1120 | 10x |
a.prop.q.g2 <- dat$param$a.prop.q.g2 |
1121 | 10x |
if (!is.null(a.prop.q.g2) && |
1122 | 10x |
!(length(a.prop.q.g2) == 1 || length(a.prop.q.g2 == nsteps))) { |
1123 | ! |
stop("Length of a.prop.q.g2 must be 1 or the value of nsteps") |
1124 |
} |
|
1125 |
|
|
1126 |
# Process ----------------------------------------------------------------- |
|
1127 | 10x |
nArrivals <- nArrivals.g2 <- 0 |
1128 |
|
|
1129 | 10x |
if (groups == 1) { |
1130 | 10x |
if (a.rand == TRUE) { |
1131 | 10x |
nArrivals <- sum(rbinom(nOld, 1, a.rate)) |
1132 |
} |
|
1133 | 10x |
if (a.rand == FALSE) { |
1134 | ! |
nArrivals <- round(nOld * a.rate) |
1135 |
} |
|
1136 |
} |
|
1137 | 10x |
if (groups == 2) { |
1138 | ! |
nOldG2 <- dat$epi$num.g2[at - 1] |
1139 | ! |
if (a.rand == TRUE) { |
1140 | ! |
if (is.na(a.rate.g2)) { |
1141 | ! |
nArrivals <- sum(rbinom(nOld, 1, a.rate)) |
1142 | ! |
nArrivals.g2 <- sum(rbinom(nOld, 1, a.rate)) |
1143 |
} else { |
|
1144 | ! |
nArrivals <- sum(rbinom(nOld, 1, a.rate)) |
1145 | ! |
nArrivals.g2 <- sum(rbinom(nOldG2, 1, a.rate.g2)) |
1146 |
} |
|
1147 |
} |
|
1148 | ! |
if (a.rand == FALSE) { |
1149 | ! |
if (is.na(a.rate.g2)) { |
1150 | ! |
nArrivals <- round(nOld * a.rate) |
1151 | ! |
nArrivals.g2 <- round(nOld * a.rate) |
1152 |
} else { |
|
1153 | ! |
nArrivals <- round(nOld * a.rate) |
1154 | ! |
nArrivals.g2 <- round(nOldG2 * a.rate.g2) |
1155 |
} |
|
1156 |
} |
|
1157 |
} |
|
1158 |
|
|
1159 |
|
|
1160 |
## Set attributes |
|
1161 | 10x |
totArrivals <- 0 |
1162 | 10x |
totArrivals.g2 <- 0 |
1163 |
|
|
1164 |
# partition arrivals into compartments |
|
1165 | 10x |
if (length(a.prop.e) > 1) { |
1166 | ! |
nArrivals.e <- round(nArrivals*(a.prop.e[at])) |
1167 | ! |
totArrivals <- totArrivals + nArrivals.e |
1168 | ! |
if (!is.null(a.prop.e.g2)) { |
1169 | ! |
nArrivals.e.g2 <- round(nArrivals.g2*(a.prop.e.g2[at])) |
1170 | ! |
totArrivals.g2 <- totArrivals.g2 + nArrivals.e.g2 |
1171 |
} else { |
|
1172 | ! |
nArrivals.e.g2 <- round(nArrivals.g2*(a.prop.e.g2[at])) |
1173 | ! |
totArrivals.g2 <- totArrivals.g2 + nArrivals.e.g2 |
1174 |
} |
|
1175 |
} else { |
|
1176 | 10x |
nArrivals.e <- round(nArrivals*a.prop.e) |
1177 | 10x |
totArrivals <- totArrivals + nArrivals.e |
1178 | 10x |
if (!is.null(a.prop.e.g2)) { |
1179 | ! |
nArrivals.e.g2 <- round(nArrivals.g2*(a.prop.e.g2)) |
1180 | ! |
totArrivals.g2 <- totArrivals.g2 + nArrivals.e.g2 |
1181 |
} else { |
|
1182 | 10x |
nArrivals.e.g2 <- round(nArrivals.g2*(a.prop.e.g2)) |
1183 | 10x |
totArrivals.g2 <- totArrivals.g2 + nArrivals.e.g2 |
1184 |
} |
|
1185 |
} |
|
1186 |
|
|
1187 | 10x |
if (length(a.prop.i) > 1) { |
1188 | ! |
nArrivals.i <- round(nArrivals*(a.prop.i[at])) |
1189 | ! |
totArrivals <- totArrivals + nArrivals.i |
1190 | ! |
if (!is.null(a.prop.i.g2)) { |
1191 | ! |
nArrivals.i.g2 <- round(nArrivals.g2*(a.prop.i.g2[at])) |
1192 | ! |
totArrivals.g2 <- totArrivals.g2 + nArrivals.i.g2 |
1193 |
} else { |
|
1194 | ! |
nArrivals.i.g2 <- round(nArrivals.g2*(a.prop.i.g2[at])) |
1195 | ! |
totArrivals.g2 <- totArrivals.g2 + nArrivals.i.g2 |
1196 |
} |
|
1197 |
} else { |
|
1198 | 10x |
nArrivals.i <- round(nArrivals*a.prop.i) |
1199 | 10x |
totArrivals <- totArrivals + nArrivals.i |
1200 | 10x |
if (!is.null(a.prop.i.g2)) { |
1201 | ! |
nArrivals.i.g2 <- round(nArrivals.g2*(a.prop.i.g2)) |
1202 | ! |
totArrivals.g2 <- totArrivals.g2 + nArrivals.i.g2 |
1203 |
} else { |
|
1204 | 10x |
nArrivals.i.g2 <- round(nArrivals.g2*(a.prop.i.g2)) |
1205 | 10x |
totArrivals.g2 <- totArrivals.g2 + nArrivals.i.g2 |
1206 |
} |
|
1207 |
} |
|
1208 |
|
|
1209 | 10x |
if (length(a.prop.q) > 1) { |
1210 | ! |
nArrivals.q <- round(nArrivals*(a.prop.q[at])) |
1211 | ! |
totArrivals <- totArrivals + nArrivals.q |
1212 | ! |
if (!is.null(a.prop.q.g2)) { |
1213 | ! |
nArrivals.q.g2 <- round(nArrivals.g2*(a.prop.q.g2[at])) |
1214 | ! |
totArrivals.g2 <- totArrivals.g2 + nArrivals.q.g2 |
1215 |
} else { |
|
1216 | ! |
nArrivals.q.g2 <- round(nArrivals.g2*(a.prop.q.g2[at])) |
1217 | ! |
totArrivals.g2 <- totArrivals.g2 + nArrivals.q.g2 |
1218 |
} |
|
1219 |
} else { |
|
1220 | 10x |
nArrivals.q <- round(nArrivals*a.prop.q) |
1221 | 10x |
totArrivals <- totArrivals + nArrivals.q |
1222 | 10x |
if (!is.null(a.prop.q.g2)) { |
1223 | ! |
nArrivals.q.g2 <- round(nArrivals.g2*(a.prop.q.g2)) |
1224 | ! |
totArrivals.g2 <- totArrivals.g2 + nArrivals.q.g2 |
1225 |
} else { |
|
1226 | 10x |
nArrivals.q.g2 <- round(nArrivals.g2*(a.prop.q.g2)) |
1227 | 10x |
totArrivals.g2 <- totArrivals.g2 + nArrivals.q.g2 |
1228 |
} |
|
1229 |
} |
|
1230 |
|
|
1231 |
# debug |
|
1232 |
# print("totArrivals:") |
|
1233 |
# print(totArrivals) |
|
1234 |
# print("totArrivals.g2:") |
|
1235 |
# print(totArrivals.g2) |
|
1236 |
# print("----") |
|
1237 |
# |
|
1238 |
# group 1 |
|
1239 | 10x |
dat$attr$active <- c(dat$attr$active, rep(1, totArrivals)) |
1240 | 10x |
dat$attr$group <- c(dat$attr$group, rep(1, totArrivals)) |
1241 | 10x |
dat$attr$status <- c(dat$attr$status, |
1242 | 10x |
rep("e", nArrivals.e), |
1243 | 10x |
rep("i", nArrivals.i), |
1244 | 10x |
rep("q", nArrivals.q), |
1245 | 10x |
rep("s", totArrivals - nArrivals.e - nArrivals.i - nArrivals.q)) |
1246 | 10x |
dat$attr$expTime <- c(dat$attr$expTime, rep(NA, totArrivals)) |
1247 | 10x |
dat$attr$infTime <- c(dat$attr$infTime, rep(NA, totArrivals)) |
1248 | 10x |
dat$attr$quarTime <- c(dat$attr$quarTime, rep(NA, totArrivals)) |
1249 | 10x |
dat$attr$hospTime <- c(dat$attr$ihospTime, rep(NA, totArrivals)) |
1250 | 10x |
dat$attr$recovTime <- c(dat$attr$recovTime, rep(NA, totArrivals)) |
1251 | 10x |
dat$attr$fatTime <- c(dat$attr$fatTime, rep(NA, totArrivals)) |
1252 |
|
|
1253 |
# group 2 |
|
1254 | 10x |
if (length(totArrivals.g2) > 0) { |
1255 | ! |
dat$attr$active <- c(dat$attr$active, rep(1, totArrivals.g2)) |
1256 | ! |
dat$attr$group <- c(dat$attr$group, rep(2, totArrivals.g2)) |
1257 | ! |
dat$attr$status <- c(dat$attr$status, |
1258 | ! |
rep("e", nArrivals.e.g2), |
1259 | ! |
rep("i", nArrivals.i.g2), |
1260 | ! |
rep("q", nArrivals.q.g2), |
1261 | ! |
rep("s", totArrivals.g2 - nArrivals.e.g2 - |
1262 | ! |
nArrivals.i.g2 - nArrivals.q.g2)) |
1263 | ! |
dat$attr$expTime <- c(dat$attr$expTime, rep(NA, totArrivals.g2)) |
1264 | ! |
dat$attr$infTime <- c(dat$attr$infTime, rep(NA, totArrivals.g2)) |
1265 | ! |
dat$attr$quarTime <- c(dat$attr$quarTime, rep(NA, totArrivals.g2)) |
1266 | ! |
dat$attr$hospTime <- c(dat$attr$ihospTime, rep(NA, totArrivals.g2)) |
1267 | ! |
dat$attr$recovTime <- c(dat$attr$recovTime, rep(NA, totArrivals.g2)) |
1268 | ! |
dat$attr$fatTime <- c(dat$attr$fatTime, rep(NA, totArrivals.g2)) |
1269 |
} |
|
1270 |
|
|
1271 |
# Output ------------------------------------------------------------------ |
|
1272 | 10x |
if (at == 2) { |
1273 | 10x |
dat$epi$a.flow <- c(0, totArrivals) |
1274 | 10x |
dat$epi$a.e.flow <- c(0, nArrivals.e) |
1275 | 10x |
dat$epi$a.i.flow <- c(0, nArrivals.i) |
1276 | 10x |
dat$epi$a.q.flow <- c(0, nArrivals.q) |
1277 |
} else { |
|
1278 | ! |
dat$epi$a.flow[at] <- totArrivals |
1279 | ! |
dat$epi$a.e.flow[at] <- c(0, nArrivals.e) |
1280 | ! |
dat$epi$a.i.flow[at] <- c(0, nArrivals.i) |
1281 | ! |
dat$epi$a.q.flow[at] <- c(0, nArrivals.q) |
1282 |
} |
|
1283 | 10x |
if (length(totArrivals.g2) > 0 && dat$param$groups == 2) { |
1284 | ! |
if (at == 2) { |
1285 | ! |
dat$epi$a.flow.g2 <- c(0, totArrivals.g2) |
1286 | ! |
dat$epi$a.e.flow.g2 <- c(0, nArrivals.e.g2) |
1287 | ! |
dat$epi$a.i.flow.g2 <- c(0, nArrivals.i.g2) |
1288 | ! |
dat$epi$a.q.flow.g2 <- c(0, nArrivals.q.g2) |
1289 |
} else { |
|
1290 | ! |
dat$epi$a.flow.g2[at] <- totArrivals.g2 |
1291 | ! |
dat$epi$a.e.flow.g2[at] <- c(0, nArrivals.e.g2) |
1292 | ! |
dat$epi$a.i.flow.g2[at] <- c(0, nArrivals.i.g2) |
1293 | ! |
dat$epi$a.q.flow.g2[at] <- c(0, nArrivals.q.g2) |
1294 |
} |
|
1295 |
} |
|
1296 |
|
|
1297 | 10x |
return(dat) |
1298 |
} |
|
1299 | ||
1300 |
departures.seiqhrf.icm <- function(dat, at, seed = NULL) { |
|
1301 |
|
|
1302 | 10x |
if(!is.null(seed)) set.seed(seed) |
1303 |
|
|
1304 |
# Conditions -------------------------------------------------------------- |
|
1305 | 10x |
if (dat$param$vital == FALSE) { |
1306 | ! |
return(dat) |
1307 |
} |
|
1308 |
|
|
1309 |
|
|
1310 |
# Variables --------------------------------------------------------------- |
|
1311 | 10x |
groups <- dat$param$groups |
1312 | 10x |
group <- dat$attr$group |
1313 |
|
|
1314 |
|
|
1315 |
# Susceptible departures ------------------------------------------------------ |
|
1316 | 10x |
nDepartures <- nDeparturesG2 <- 0 |
1317 | 10x |
idsElig <- which(dat$attr$active == 1 & dat$attr$status == "s") |
1318 | 10x |
nElig <- length(idsElig) |
1319 | 10x |
if (nElig > 0) { |
1320 |
|
|
1321 | 10x |
gElig <- group[idsElig] |
1322 | 10x |
rates <- c(dat$param$ds.rate, dat$param$ds.rate.g2) |
1323 | 10x |
ratesElig <- rates[gElig] |
1324 |
|
|
1325 | 10x |
if (dat$control$d.rand == TRUE) { |
1326 | 10x |
vecDepartures <- which(rbinom(nElig, 1, ratesElig) == 1) |
1327 | 10x |
if (length(vecDepartures) > 0) { |
1328 | 1x |
idsDpt <- idsElig[vecDepartures] |
1329 | 1x |
nDepartures <- sum(group[idsDpt] == 1) |
1330 | 1x |
nDeparturesG2 <- sum(group[idsDpt] == 2) |
1331 | 1x |
dat$attr$active[idsDpt] <- 0 |
1332 |
} |
|
1333 |
} else { |
|
1334 | ! |
nDepartures <- min(round(sum(ratesElig[gElig == 1])), sum(gElig == 1)) |
1335 | ! |
dat$attr$active[ssample(idsElig[gElig == 1], nDepartures)] <- 0 |
1336 | ! |
if (groups == 2) { |
1337 | ! |
nDeparturesG2 <- min(round(sum(ratesElig[gElig == 2])), sum(gElig == 2)) |
1338 | ! |
dat$attr$active[ssample(idsElig[gElig == 2], nDeparturesG2)] <- 0 |
1339 |
} |
|
1340 |
} |
|
1341 |
} |
|
1342 |
|
|
1343 | 10x |
if (at == 2) { |
1344 | 10x |
dat$epi$ds.flow <- c(0, nDepartures) |
1345 | 10x |
if (groups == 2) { |
1346 | ! |
dat$epi$ds.flow.g2 <- c(0, nDeparturesG2) |
1347 |
} |
|
1348 |
} else { |
|
1349 | ! |
dat$epi$ds.flow[at] <- nDepartures |
1350 | ! |
if (groups == 2) { |
1351 | ! |
dat$epi$ds.flow.g2[at] <- nDeparturesG2 |
1352 |
} |
|
1353 |
} |
|
1354 |
|
|
1355 |
# Exposed Departures --------------------------------------------------------- |
|
1356 | 10x |
nDepartures <- nDeparturesG2 <- 0 |
1357 | 10x |
idsElig <- which(dat$attr$active == 1 & dat$attr$status == "e") |
1358 | 10x |
nElig <- length(idsElig) |
1359 | 10x |
if (nElig > 0) { |
1360 |
|
|
1361 | 10x |
gElig <- group[idsElig] |
1362 | 10x |
rates <- c(dat$param$de.rate, dat$param$de.rate.g2) |
1363 | 10x |
ratesElig <- rates[gElig] |
1364 |
|
|
1365 | 10x |
if (dat$control$d.rand == TRUE) { |
1366 | 10x |
vecDepartures <- which(rbinom(nElig, 1, ratesElig) == 1) |
1367 | 10x |
if (length(vecDepartures) > 0) { |
1368 | ! |
idsDpt <- idsElig[vecDepartures] |
1369 | ! |
nDepartures <- sum(group[idsDpt] == 1) |
1370 | ! |
nDeparturesG2 <- sum(group[idsDpt] == 2) |
1371 | ! |
dat$attr$active[idsDpt] <- 0 |
1372 |
} |
|
1373 |
} else { |
|
1374 | ! |
nDepartures <- min(round(sum(ratesElig[gElig == 1])), sum(gElig == 1)) |
1375 | ! |
dat$attr$active[ssample(idsElig[gElig == 1], nDepartures)] <- 0 |
1376 | ! |
if (groups == 2) { |
1377 | ! |
nDeparturesG2 <- min(round(sum(ratesElig[gElig == 2])), sum(gElig == 2)) |
1378 | ! |
dat$attr$active[ssample(idsElig[gElig == 2], nDeparturesG2)] <- 0 |
1379 |
} |
|
1380 |
} |
|
1381 |
} |
|
1382 |
|
|
1383 | 10x |
if (at == 2) { |
1384 | 10x |
dat$epi$de.flow <- c(0, nDepartures) |
1385 | 10x |
if (groups == 2) { |
1386 | ! |
dat$epi$de.flow.g2 <- c(0, nDeparturesG2) |
1387 |
} |
|
1388 |
} else { |
|
1389 | ! |
dat$epi$de.flow[at] <- nDepartures |
1390 | ! |
if (groups == 2) { |
1391 | ! |
dat$epi$de.flow.g2[at] <- nDeparturesG2 |
1392 |
} |
|
1393 |
} |
|
1394 |
|
|
1395 |
|
|
1396 |
# Infected Departures --------------------------------------------------------- |
|
1397 | 10x |
nDepartures <- nDeparturesG2 <- 0 |
1398 | 10x |
idsElig <- which(dat$attr$active == 1 & dat$attr$status == "i") |
1399 | 10x |
nElig <- length(idsElig) |
1400 | 10x |
if (nElig > 0) { |
1401 |
|
|
1402 | 10x |
gElig <- group[idsElig] |
1403 | 10x |
rates <- c(dat$param$di.rate, dat$param$di.rate.g2) |
1404 | 10x |
ratesElig <- rates[gElig] |
1405 |
|
|
1406 | 10x |
if (dat$control$d.rand == TRUE) { |
1407 | 10x |
vecDepartures <- which(rbinom(nElig, 1, ratesElig) == 1) |
1408 | 10x |
if (length(vecDepartures) > 0) { |
1409 | ! |
idsDpt <- idsElig[vecDepartures] |
1410 | ! |
nDepartures <- sum(group[idsDpt] == 1) |
1411 | ! |
nDeparturesG2 <- sum(group[idsDpt] == 2) |
1412 | ! |
dat$attr$active[idsDpt] <- 0 |
1413 |
} |
|
1414 |
} else { |
|
1415 | ! |
nDepartures <- min(round(sum(ratesElig[gElig == 1])), sum(gElig == 1)) |
1416 | ! |
dat$attr$active[ssample(idsElig[gElig == 1], nDepartures)] <- 0 |
1417 | ! |
if (groups == 2) { |
1418 | ! |
nDeparturesG2 <- min(round(sum(ratesElig[gElig == 2])), sum(gElig == 2)) |
1419 | ! |
dat$attr$active[ssample(idsElig[gElig == 2], nDeparturesG2)] <- 0 |
1420 |
} |
|
1421 |
} |
|
1422 |
} |
|
1423 |
|
|
1424 | 10x |
if (at == 2) { |
1425 | 10x |
dat$epi$di.flow <- c(0, nDepartures) |
1426 | 10x |
if (groups == 2) { |
1427 | ! |
dat$epi$di.flow.g2 <- c(0, nDeparturesG2) |
1428 |
} |
|
1429 |
} else { |
|
1430 | ! |
dat$epi$di.flow[at] <- nDepartures |
1431 | ! |
if (groups == 2) { |
1432 | ! |
dat$epi$di.flow.g2[at] <- nDeparturesG2 |
1433 |
} |
|
1434 |
} |
|
1435 |
|
|
1436 |
# Quarantined Departures --------------------------------------------------------- |
|
1437 | 10x |
nDepartures <- nDeparturesG2 <- 0 |
1438 | 10x |
idsElig <- which(dat$attr$active == 1 & dat$attr$status == "q") |
1439 | 10x |
nElig <- length(idsElig) |
1440 | 10x |
if (nElig > 0) { |
1441 |
|
|
1442 | ! |
gElig <- group[idsElig] |
1443 | ! |
rates <- c(dat$param$dq.rate, dat$param$dq.rate.g2) |
1444 | ! |
ratesElig <- rates[gElig] |
1445 |
|
|
1446 | ! |
if (dat$control$d.rand == TRUE) { |
1447 | ! |
vecDepartures <- which(rbinom(nElig, 1, ratesElig) == 1) |
1448 | ! |
if (length(vecDepartures) > 0) { |
1449 | ! |
idsDpt <- idsElig[vecDepartures] |
1450 | ! |
nDepartures <- sum(group[idsDpt] == 1) |
1451 | ! |
nDeparturesG2 <- sum(group[idsDpt] == 2) |
1452 | ! |
dat$attr$active[idsDpt] <- 0 |
1453 |
} |
|
1454 |
} else { |
|
1455 | ! |
nDepartures <- min(round(sum(ratesElig[gElig == 1])), sum(gElig == 1)) |
1456 | ! |
dat$attr$active[ssample(idsElig[gElig == 1], nDepartures)] <- 0 |
1457 | ! |
if (groups == 2) { |
1458 | ! |
nDeparturesG2 <- min(round(sum(ratesElig[gElig == 2])), sum(gElig == 2)) |
1459 | ! |
dat$attr$active[ssample(idsElig[gElig == 2], nDeparturesG2)] <- 0 |
1460 |
} |
|
1461 |
} |
|
1462 |
} |
|
1463 |
|
|
1464 | 10x |
if (at == 2) { |
1465 | 10x |
dat$epi$dq.flow <- c(0, nDepartures) |
1466 | 10x |
if (groups == 2) { |
1467 | ! |
dat$epi$dq.flow.g2 <- c(0, nDeparturesG2) |
1468 |
} |
|
1469 |
} else { |
|
1470 | ! |
dat$epi$dq.flow[at] <- nDepartures |
1471 | ! |
if (groups == 2) { |
1472 | ! |
dat$epi$dq.flow.g2[at] <- nDeparturesG2 |
1473 |
} |
|
1474 |
} |
|
1475 |
|
|
1476 |
# Hospitalised Departures --------------------------------------------------------- |
|
1477 | 10x |
nDepartures <- nDeparturesG2 <- 0 |
1478 | 10x |
idsElig <- which(dat$attr$active == 1 & dat$attr$status == "h") |
1479 | 10x |
nElig <- length(idsElig) |
1480 | 10x |
if (nElig > 0) { |
1481 |
|
|
1482 | ! |
gElig <- group[idsElig] |
1483 | ! |
rates <- c(dat$param$dh.rate, dat$param$dh.rate.g2) |
1484 | ! |
ratesElig <- rates[gElig] |
1485 |
|
|
1486 | ! |
if (dat$control$d.rand == TRUE) { |
1487 | ! |
vecDepartures <- which(rbinom(nElig, 1, ratesElig) == 1) |
1488 | ! |
if (length(vecDepartures) > 0) { |
1489 | ! |
idsDpt <- idsElig[vecDepartures] |
1490 | ! |
nDepartures <- sum(group[idsDpt] == 1) |
1491 | ! |
nDeparturesG2 <- sum(group[idsDpt] == 2) |
1492 | ! |
dat$attr$active[idsDpt] <- 0 |
1493 |
} |
|
1494 |
} else { |
|
1495 | ! |
nDepartures <- min(round(sum(ratesElig[gElig == 1])), sum(gElig == 1)) |
1496 | ! |
dat$attr$active[ssample(idsElig[gElig == 1], nDepartures)] <- 0 |
1497 | ! |
if (groups == 2) { |
1498 | ! |
nDeparturesG2 <- min(round(sum(ratesElig[gElig == 2])), sum(gElig == 2)) |
1499 | ! |
dat$attr$active[ssample(idsElig[gElig == 2], nDeparturesG2)] <- 0 |
1500 |
} |
|
1501 |
} |
|
1502 |
} |
|
1503 |
|
|
1504 | 10x |
if (at == 2) { |
1505 | 10x |
dat$epi$dh.flow <- c(0, nDepartures) |
1506 | 10x |
if (groups == 2) { |
1507 | ! |
dat$epi$dh.flow.g2 <- c(0, nDeparturesG2) |
1508 |
} |
|
1509 |
} else { |
|
1510 | ! |
dat$epi$dh.flow[at] <- nDepartures |
1511 | ! |
if (groups == 2) { |
1512 | ! |
dat$epi$dh.flow.g2[at] <- nDeparturesG2 |
1513 |
} |
|
1514 |
} |
|
1515 |
|
|
1516 |
|
|
1517 |
# Recovered Departures -------------------------------------------------------- |
|
1518 | 10x |
nDepartures <- nDeparturesG2 <- 0 |
1519 | 10x |
idsElig <- which(dat$attr$active == 1 & dat$attr$status == "r") |
1520 | 10x |
nElig <- length(idsElig) |
1521 | 10x |
if (nElig > 0) { |
1522 |
|
|
1523 | ! |
gElig <- group[idsElig] |
1524 | ! |
rates <- c(dat$param$dr.rate, dat$param$dr.rate.g2) |
1525 | ! |
ratesElig <- rates[gElig] |
1526 |
|
|
1527 | ! |
if (dat$control$d.rand == TRUE) { |
1528 | ! |
vecDepartures <- which(rbinom(nElig, 1, ratesElig) == 1) |
1529 | ! |
if (length(vecDepartures) > 0) { |
1530 | ! |
idsDpt <- idsElig[vecDepartures] |
1531 | ! |
nDepartures <- sum(group[idsDpt] == 1) |
1532 | ! |
nDeparturesG2 <- sum(group[idsDpt] == 2) |
1533 | ! |
dat$attr$active[idsDpt] <- 0 |
1534 |
} |
|
1535 |
} else { |
|
1536 | ! |
nDepartures <- min(round(sum(ratesElig[gElig == 1])), sum(gElig == 1)) |
1537 | ! |
dat$attr$active[ssample(idsElig[gElig == 1], nDepartures)] <- 0 |
1538 | ! |
if (groups == 2) { |
1539 | ! |
nDeparturesG2 <- min(round(sum(ratesElig[gElig == 2])), sum(gElig == 2)) |
1540 | ! |
dat$attr$active[ssample(idsElig[gElig == 2], nDeparturesG2)] <- 0 |
1541 |
} |
|
1542 |
} |
|
1543 |
} |
|
1544 |
|
|
1545 | 10x |
if (at == 2) { |
1546 | 10x |
dat$epi$dr.flow <- c(0, nDepartures) |
1547 | 10x |
if (groups == 2) { |
1548 | ! |
dat$epi$dr.flow.g2 <- c(0, nDeparturesG2) |
1549 |
} |
|
1550 |
} else { |
|
1551 | ! |
dat$epi$dr.flow[at] <- nDepartures |
1552 | ! |
if (groups == 2) { |
1553 | ! |
dat$epi$dr.flow.g2[at] <- nDeparturesG2 |
1554 |
} |
|
1555 |
} |
|
1556 |
|
|
1557 | 10x |
return(dat) |
1558 |
} |
1 |
## internal function |
|
2 | ||
3 |
get_del <- function(dat, lab, acts, letgo = FALSE, seed = NULL){ |
|
4 |
|
|
5 |
#if(!is.null(seed)) set.seed(seed) |
|
6 |
|
|
7 |
## Edgelist |
|
8 | 36x |
p1 <- EpiModel::ssample(which(dat$attr$active == 1 & dat$attr$status != "f"), acts, replace = TRUE) ## Bootstrap index, allow repetitive elements |
9 | 36x |
p2 <- EpiModel::ssample(which(dat$attr$active == 1 & dat$attr$status != "f"), acts, replace = TRUE) |
10 |
|
|
11 |
#### keep samplling until all elements in p1 and p2 differ |
|
12 | 36x |
del <- NULL |
13 | 36x |
if (length(p1) > 0 & length(p2) > 0) { |
14 | 36x |
del <- data.frame(p1, p2) |
15 | 36x |
while (any(del$p1 == del$p2)) { |
16 | 31x |
del$p2 <- ifelse(del$p1 == del$p2, |
17 | 31x |
EpiModel::ssample(which(dat$attr$active == 1 & dat$attr$status != "f"), 1), del$p2) |
18 |
} |
|
19 |
}else{ |
|
20 | ! |
if(letgo) return(NULL) |
21 |
} |
|
22 |
|
|
23 |
## Discordant edgelist (del) |
|
24 | 36x |
del$p1.stat <- dat$attr$status[del$p1] |
25 | 36x |
del$p2.stat <- dat$attr$status[del$p2] |
26 | 36x |
serodis <- (del$p1.stat == lab[1] & del$p2.stat == lab[2]) | |
27 | 36x |
(del$p1.stat == lab[2] & del$p2.stat == lab[1]) |
28 |
|
|
29 | 36x |
return(del[serodis, ]) |
30 |
|
|
31 |
} |
|
32 | ||
33 | ||
34 |
#' Infection function |
|
35 |
#' |
|
36 |
#' Function to implement infection processes for SEIQHRF model |
|
37 |
#' |
|
38 |
#' @param dat Merged input parameters. |
|
39 |
#' @param at Step number |
|
40 |
#' @param seed random seed |
|
41 |
#' |
|
42 |
#' @return Updated dat |
|
43 |
#' |
|
44 |
#' @importFrom EpiModel ssample |
|
45 |
#' @importFrom stats rbinom |
|
46 |
#' @importFrom dplyr between |
|
47 |
#' @export |
|
48 |
infection.FUN <- function(dat, at, seed = NULL){ |
|
49 |
|
|
50 | 10x |
if(!is.null(seed)) set.seed(seed) |
51 |
|
|
52 | 12x |
type <- dat$control$type |
53 | 12x |
nsteps <- dat$control$nsteps |
54 |
|
|
55 |
## Extract param |
|
56 | 12x |
act.rate.i <- dat$param$act.rate.i |
57 | 12x |
inf.prob.i <- dat$param$inf.prob.i |
58 | 12x |
if(type %in% c("SEIQHR", "SEIQHRF")){ |
59 | 12x |
quar.rate <- dat$param$quar.rate |
60 | 12x |
disch.rate <- dat$param$disch.rate |
61 | 12x |
act.rate.e <- dat$param$act.rate.e |
62 | 12x |
act.rate.q <- dat$param$act.rate.q |
63 | 12x |
inf.prob.e <- dat$param$inf.prob.e |
64 | 12x |
inf.prob.q <- dat$param$inf.prob.q |
65 |
} |
|
66 |
|
|
67 |
# Transmission from infected (a) |
|
68 |
## Expected acts |
|
69 | 12x |
acts <- round(ifelse(length(act.rate.i) > 1, act.rate.i[at - 1], act.rate.i) * |
70 | 12x |
dat$epi$num[at - 1] / 2 ) |
71 |
|
|
72 |
## Edgelist |
|
73 | 12x |
del <- get_del(dat, c("s", "i"), acts) |
74 |
|
|
75 |
## Transmission on edgelist |
|
76 | 12x |
if (nrow(del) > 0) { |
77 |
|
|
78 | 12x |
del$tprob <- ifelse(length(inf.prob.i) > 1, inf.prob.i[at], inf.prob.i) |
79 |
|
|
80 | 12x |
if (!is.null(dat$param$inter.eff.i) && |
81 | 12x |
dplyr::between(at, dat$param$inter.start.i, dat$param$inter.stop.i)) { |
82 | ! |
del$tprob <- del$tprob * (1 - dat$param$inter.eff.i) |
83 |
} |
|
84 |
|
|
85 | 12x |
del$trans <- stats::rbinom(nrow(del), 1, del$tprob) |
86 | 12x |
del <- del[del$trans == TRUE, ] |
87 |
|
|
88 | 12x |
if (nrow(del) > 0) { |
89 | 10x |
newIds <- unique(ifelse(del$p1.stat == "s", del$p1, del$p2)) ## estract individual ID for "s" |
90 | 10x |
nExp.i <- length(newIds) |
91 | 10x |
dat$attr$status[newIds] <- "e" |
92 | 10x |
dat$attr$expTime[newIds] <- at |
93 |
} else { |
|
94 | 2x |
nExp.i <- 0 |
95 |
} |
|
96 |
|
|
97 |
} else { |
|
98 | ! |
nExp.i <- 0 |
99 |
} |
|
100 |
|
|
101 | 12x |
if (type == "SEIQHRF"){ |
102 |
|
|
103 |
# Transmission from exposed |
|
104 |
## Expected acts |
|
105 | 12x |
acts <- round(ifelse(length(act.rate.e) > 1, act.rate.e[at - 1], act.rate.e) * |
106 | 12x |
dat$epi$num[at - 1] / 2 ) |
107 |
|
|
108 |
## Edgelist: s, e |
|
109 | 12x |
del <- get_del(dat, c("s", "e"), acts, letgo = TRUE) |
110 | 12x |
if(is.null(del)){ |
111 | ! |
nExp.e <- 0 |
112 |
}else{ |
|
113 |
## Transmission on edgelist |
|
114 | 12x |
if (nrow(del) > 0) { |
115 |
|
|
116 | 10x |
del$tprob <- ifelse(length(inf.prob.e) > 1, inf.prob.e[at], inf.prob.e) |
117 |
|
|
118 | 10x |
if (!is.null(dat$param$inter.eff.e) && |
119 | 10x |
dplyr::between(at, dat$param$inter.start.e, dat$param$inter.stop.e)) { |
120 | ! |
del$tprob <- del$tprob * (1 - dat$param$inter.eff.e) |
121 |
} |
|
122 |
|
|
123 | 10x |
del$trans <- stats::rbinom(nrow(del), 1, del$tprob) |
124 | 10x |
del <- del[del$trans == TRUE, ] |
125 |
|
|
126 | 10x |
if (nrow(del) > 0) { |
127 | 4x |
newIds <- unique(ifelse(del$p1.stat == "s", del$p1, del$p2)) |
128 | 4x |
nExp.e <- length(newIds) |
129 | 4x |
dat$attr$status[newIds] <- "e" |
130 | 4x |
dat$attr$expTime[newIds] <- at |
131 |
} else { |
|
132 | 6x |
nExp.e <- 0 |
133 |
} |
|
134 |
|
|
135 |
} else { |
|
136 | 2x |
nExp.e <- 0 |
137 |
} |
|
138 |
} |
|
139 |
|
|
140 |
# Transmission from quarantined |
|
141 |
## Expected acts |
|
142 | 12x |
acts <- round(ifelse(length(act.rate.q) > 1, act.rate.q[at - 1], act.rate.q) * |
143 | 12x |
dat$epi$num[at - 1] / 2 ) |
144 |
|
|
145 |
## Edgelist:s, q |
|
146 | 12x |
del <- get_del(dat, c("s", "q"), acts, letgo = TRUE, seed) |
147 | 12x |
if(is.null(del)){ |
148 | ! |
nExp.q <- 0 |
149 |
}else{ |
|
150 | 12x |
if (nrow(del) > 0) { |
151 |
|
|
152 | ! |
del$tprob <- ifelse(length(inf.prob.q) > 1, inf.prob.q[at], inf.prob.q) |
153 |
|
|
154 | ! |
if (!is.null(dat$param$inter.eff.q) && |
155 | ! |
dplyr::between(at, dat$param$inter.start.q, dat$param$inter.stop.q)) { |
156 | ! |
del$tprob <- del$tprob * (1 - dat$param$inter.eff.q) |
157 |
} |
|
158 |
|
|
159 | ! |
del$trans <- stats::rbinom(nrow(del), 1, del$tprob) |
160 | ! |
del <- del[del$trans == TRUE, ] |
161 | ! |
if (nrow(del) > 0) { |
162 | ! |
newIds <- unique(ifelse(del$p1.stat == "s", del$p1, del$p2)) |
163 | ! |
nExp.q <- length(newIds) |
164 | ! |
dat$attr$status[newIds] <- "e" |
165 | ! |
dat$attr$expTime[newIds] <- at |
166 |
} else { |
|
167 | ! |
nExp.q <- 0 |
168 |
} |
|
169 |
} else { |
|
170 | 12x |
nExp.q <- 0 |
171 |
} |
|
172 |
} |
|
173 |
|
|
174 |
} |
|
175 |
|
|
176 |
## Output |
|
177 | 12x |
tmp <- ifelse(type %in% c("SEIQHR", "SEIQHRF"), nExp.i + nExp.q, nExp.i) |
178 | 12x |
if (at == 2) { |
179 | 12x |
dat$epi$se.flow <- c(0, tmp) |
180 |
} else { |
|
181 | ! |
dat$epi$se.flow[at] <- tmp |
182 |
} |
|
183 |
|
|
184 | 12x |
dat |
185 |
} |
1 |
#' Merge icm |
|
2 |
#' |
|
3 |
#' Function to merge icms |
|
4 |
#' |
|
5 |
#' @param x First to merge |
|
6 |
#' @param y Second to merge |
|
7 |
#' @param ... other parameters |
|
8 |
#' |
|
9 |
#' @return Combined |
|
10 |
#' |
|
11 |
#' @importFrom doParallel registerDoParallel |
|
12 |
#' |
|
13 |
merge.seiqhrf.icm <- function(x, y, ...) { |
|
14 | ||
15 |
## Check structure |
|
16 | 14x |
if (length(x) != length(y) || !identical(names(x), names(y))) { |
17 | ! |
stop("x and y have different structure") |
18 |
} |
|
19 | 14x |
if (x$control$nsims > 1 & y$control$nsims > 1 & |
20 | 14x |
!identical(sapply(x, class), sapply(y, class))) { |
21 | ! |
stop("x and y have different structure") |
22 |
} |
|
23 | ||
24 |
## Check params |
|
25 | 14x |
check1 <- identical(x$param, y$param) |
26 | 14x |
check2 <- identical(x$control[-which(names(x$control) %in% c("nsims", "currsim"))], |
27 | 14x |
y$control[-which(names(y$control) %in% c("nsims", "currsim"))]) |
28 | ||
29 | 14x |
if (check1 == FALSE) { |
30 | ! |
stop("x and y have different parameters") |
31 |
} |
|
32 | 14x |
if (check2 == FALSE) { |
33 | ! |
stop("x and y have different controls") |
34 |
} |
|
35 | ||
36 | ||
37 | 14x |
z <- x |
38 | 14x |
new.range <- (x$control$nsims + 1):(x$control$nsims + y$control$nsims) |
39 | ||
40 |
# Merge data |
|
41 | 14x |
for (i in 1:length(x$epi)) { |
42 | 322x |
if (x$control$nsims == 1) { |
43 | 46x |
x$epi[[i]] <- data.frame(x$epi[[i]]) |
44 |
} |
|
45 | 322x |
if (y$control$nsims == 1) { |
46 | 322x |
y$epi[[i]] <- data.frame(y$epi[[i]]) |
47 |
} |
|
48 | 322x |
z$epi[[i]] <- cbind(x$epi[[i]], y$epi[[i]]) |
49 | 322x |
names(z$epi[[i]])[new.range] <- paste0("sim", new.range) |
50 |
} |
|
51 | ||
52 | 14x |
z$control$nsims <- max(new.range) |
53 | ||
54 | 14x |
return(z) |
55 |
} |
|
56 | ||
57 |
#' SEIQHRF function |
|
58 |
#' |
|
59 |
#' Function to implement infection processes for SEIQHRF model |
|
60 |
#' |
|
61 |
#' @param init Initial state values, returned by \code{\link{init_seiqhrf}}. |
|
62 |
#' @param control Control parameters of the model, returned by \code{\link{control_seiqhrf}}. |
|
63 |
#' @param param Model parameters, returned by \code{\link{param_seiqhrf}}. |
|
64 |
#' @param seed random seed for checking consistency |
|
65 |
#' |
|
66 |
#' @return Updated dat |
|
67 |
#' |
|
68 |
#' @importFrom foreach "%dopar%" foreach |
|
69 |
#' @importFrom EpiModel verbose.icm |
|
70 |
#' @importFrom future availableCores |
|
71 |
#' @export |
|
72 |
seiqhrf <- function(init, control, param, seed = NULL) { |
|
73 | ||
74 | ||
75 | 2x |
crosscheck_seiqhrf(param, init, control) |
76 | 2x |
nsims <- control$nsims |
77 | 2x |
ncores <- ifelse(control$nsims == 1, 1, min(future::availableCores(), control$ncores)) |
78 | 2x |
control$ncores <- ncores |
79 | ||
80 | 2x |
if (ncores == 1) { |
81 | ||
82 |
# Simulation loop start |
|
83 | ! |
for (s in 1:control$nsims) { |
84 | ||
85 |
## Initialization module |
|
86 | ! |
if (!is.null(control[["initialize.FUN"]])) { |
87 | ! |
dat <- do.call(control[["initialize.FUN"]], list(param, init, control, seed)) |
88 |
} |
|
89 | ||
90 | ||
91 |
# Timestep loop |
|
92 | ! |
for (at in 2:control$nsteps) { |
93 | ||
94 |
## Infection |
|
95 | ! |
if (!is.null(control[["infection.FUN"]])) { |
96 | ! |
dat <- do.call(control[["infection.FUN"]], list(dat, at, seed)) |
97 |
} |
|
98 | ||
99 | ||
100 |
## Recovery |
|
101 | ! |
if (!is.null(control[["recovery.FUN"]])) { |
102 | ! |
dat <- do.call(control[["recovery.FUN"]], list(dat, at, seed)) |
103 |
} |
|
104 | ||
105 | ||
106 |
## Departure Module |
|
107 | ! |
if (!is.null(control[["departures.FUN"]])) { |
108 | ! |
dat <- do.call(control[["departures.FUN"]], list(dat, at, seed)) |
109 |
} |
|
110 | ||
111 | ||
112 |
## Arrival Module |
|
113 | ! |
if (!is.null(control[["arrivals.FUN"]])) { |
114 | ! |
dat <- do.call(control[["arrivals.FUN"]], list(dat, at, seed)) |
115 |
} |
|
116 | ||
117 | ||
118 |
## Outputs |
|
119 | ! |
if (!is.null(control[["get_prev.FUN"]])) { |
120 | ! |
dat <- do.call(control[["get_prev.FUN"]], list(dat, at)) |
121 |
} |
|
122 | ||
123 | ||
124 |
## Track progress |
|
125 | ! |
EpiModel::verbose.icm(dat, type = "progress", s, at) |
126 |
} |
|
127 | ||
128 |
# Set output |
|
129 | ! |
if (s == 1) { |
130 | ! |
out <- saveout_seiqhrf(dat, s) |
131 |
} else { |
|
132 | ! |
out <- saveout_seiqhrf(dat, s, out) |
133 |
} |
|
134 | ||
135 |
} # Simulation loop end |
|
136 | ||
137 |
} # end of single core execution |
|
138 | ||
139 | 2x |
if (ncores > 1) { |
140 | 2x |
doParallel::registerDoParallel(ncores) |
141 | ||
142 | 2x |
sout <- foreach(s = 1:nsims) %dopar% { |
143 | ||
144 | ! |
control$nsims <- 1 |
145 | ! |
control$currsim <- s |
146 | ||
147 |
## Initialization module |
|
148 | ! |
if (!is.null(control[["initialize.FUN"]])) { |
149 | ! |
dat <- do.call(control[["initialize.FUN"]], list(param, init, control, seed)) |
150 |
} |
|
151 | ||
152 |
# Timestep loop |
|
153 | ! |
for (at in 2:control$nsteps) { |
154 | ||
155 |
## Infection |
|
156 | ! |
if (!is.null(control[["infection.FUN"]])) { |
157 | ! |
dat <- do.call(control[["infection.FUN"]], list(dat, at, seed)) |
158 |
} |
|
159 | ||
160 | ||
161 |
## Recovery |
|
162 | ! |
if (!is.null(control[["recovery.FUN"]])) { |
163 | ! |
dat <- do.call(control[["recovery.FUN"]], list(dat, at, seed)) |
164 |
} |
|
165 | ||
166 | ||
167 |
## Departure Module |
|
168 | ! |
if (!is.null(control[["departures.FUN"]])) { |
169 | ! |
dat <- do.call(control[["departures.FUN"]], list(dat, at, seed)) |
170 |
} |
|
171 | ||
172 | ||
173 |
## Arrival Module |
|
174 | ! |
if (!is.null(control[["arrivals.FUN"]])) { |
175 | ! |
dat <- do.call(control[["arrivals.FUN"]], list(dat, at, seed)) |
176 |
} |
|
177 | ||
178 | ||
179 |
## Outputs |
|
180 | ! |
if (!is.null(control[["get_prev.FUN"]])) { |
181 | ! |
dat <- do.call(control[["get_prev.FUN"]], list(dat, at)) |
182 |
} |
|
183 | ||
184 | ||
185 |
## Track progress |
|
186 | ! |
EpiModel::verbose.icm(dat, type = "progress", s, at) |
187 |
} |
|
188 | ||
189 |
# Set output |
|
190 | ! |
out <- saveout_seiqhrf(dat, s=1) |
191 | ! |
class(out) <- "icm" |
192 | ! |
return(out) |
193 | ||
194 |
} |
|
195 | ||
196 |
# aggregate results collected from each thread |
|
197 | 2x |
collected_times <- list() |
198 | ||
199 |
# collect the times from sout then delete them |
|
200 | 2x |
for (i in 1:length(sout)) { |
201 | 16x |
collected_times[[paste0("sim", i)]] <- sout[[i]]$times$sim1 |
202 | 16x |
sout[[i]]$times <- NULL |
203 |
} |
|
204 | ||
205 |
# merge $epi structures |
|
206 | 2x |
merged.out <- sout[[1]] |
207 | 2x |
for (i in 2:length(sout)) { |
208 | 14x |
merged.out <- merge.seiqhrf.icm(merged.out, sout[[i]], param.error = FALSE) |
209 |
} |
|
210 | 2x |
out <- merged.out |
211 | ||
212 |
# add the collected timing data |
|
213 | 2x |
out$times <- collected_times |
214 | ||
215 |
} # end of parallel execution |
|
216 |
|
|
217 | 2x |
class(out) <- "seiqhrf" |
218 | ||
219 | 2x |
return(out) |
220 |
} |
1 |
#' @export |
|
2 |
print.init.seiqhrf <- function(x, ...) { |
|
3 |
|
|
4 | ! |
cat("SEIQHRF Initial Conditions") |
5 | ! |
cat("\n===========================\n") |
6 | ! |
cat("\nUser specified control parameters:") |
7 | ! |
cat("\n---------------------------\n") |
8 | ! |
for(i in x$usr.specified){ |
9 | ! |
if (class(x[[i]]) %in% c("integer", "numeric") && length(x[[i]]) > 10) { |
10 | ! |
cat(names(x)[i], "=", x[[i]][1:5], "...", fill = 80) |
11 | ! |
} else if (class(x[[i]]) == "data.frame") { |
12 | ! |
cat(names(x)[i], "= <data.frame>\n") |
13 | ! |
} else if (class(x[[i]]) == "list") { |
14 | ! |
cat(i, "= <list>\n") |
15 |
} else { |
|
16 | ! |
cat(i, "=", x[[i]], fill = 80) |
17 |
} |
|
18 |
} |
|
19 | ! |
cat("\nDefault control parameters:") |
20 | ! |
cat("\n---------------------------\n") |
21 |
|
|
22 | ! |
pToPrint <- names(x)[seq_along(x)] |
23 | ! |
pToPrint<- setdiff(pToPrint, c(x$usr.specified, "usr.specified")) |
24 | ! |
for (i in pToPrint) { |
25 | ! |
if (class(x[[i]]) %in% c("integer", "numeric") && length(x[[i]]) > 10) { |
26 | ! |
cat(names(x)[i], "=", x[[i]][1:5], "...", fill = 80) |
27 | ! |
} else if (class(x[[i]]) == "data.frame") { |
28 | ! |
cat(names(x)[i], "= <data.frame>\n") |
29 | ! |
} else if (class(x[[i]]) == "list") { |
30 | ! |
cat(i, "= <list>\n") |
31 |
} else { |
|
32 | ! |
cat(i, "=", x[[i]], fill = 80) |
33 |
} |
|
34 |
} |
|
35 |
|
|
36 | ! |
invisible() |
37 |
} |
|
38 | ||
39 | ||
40 |
#' @export |
|
41 |
print.param.seiqhrf <- function(x, ...) { |
|
42 |
|
|
43 | ! |
cat("SEIQHRF Parameters") |
44 | ! |
cat("\n===========================\n") |
45 |
|
|
46 | ! |
cat("\nUser specified control parameters:") |
47 | ! |
cat("\n---------------------------\n") |
48 | ! |
for(i in x$usr.specified){ |
49 | ! |
if (class(x[[i]]) %in% c("integer", "numeric") && length(x[[i]]) > 10) { |
50 | ! |
cat(i, "=", x[[i]][1:5], "...", fill = 80) |
51 | ! |
} else if (class(x[[i]]) == "data.frame") { |
52 | ! |
cat(i, "= <data.frame>\n") |
53 | ! |
} else if (class(x[[i]]) == "list") { |
54 | ! |
cat(i, "= <list>\n") |
55 |
} else { |
|
56 | ! |
cat(i, "=", x[[i]], fill = 80) |
57 |
} |
|
58 |
} |
|
59 |
|
|
60 | ! |
cat("\nDefault control parameters:") |
61 | ! |
cat("\n---------------------------\n") |
62 | ! |
pToPrint <- names(x)[which(!(names(x) %in% c("vital")))] |
63 | ! |
pToPrint<- setdiff(pToPrint, c(x$usr.specified, "usr.specified")) |
64 | ! |
for (i in pToPrint) { |
65 | ! |
if (class(x[[i]]) %in% c("integer", "numeric") && length(x[[i]]) > 10) { |
66 | ! |
cat(i, "=", x[[i]][1:5], "...", fill = 80) |
67 | ! |
} else if (class(x[[i]]) == "data.frame") { |
68 | ! |
cat(i, "= <data.frame>\n") |
69 | ! |
} else if (class(x[[i]]) == "list") { |
70 | ! |
cat(i, "= <list>\n") |
71 |
} else { |
|
72 | ! |
cat(i, "=", x[[i]], fill = 80) |
73 |
} |
|
74 |
} |
|
75 |
|
|
76 | ! |
invisible() |
77 |
} |
|
78 | ||
79 | ||
80 |
#' @export |
|
81 |
print.control.seiqhrf <- function(x, ...) { |
|
82 |
|
|
83 | ! |
cat("SEIQHRF Control Settings") |
84 | ! |
cat("\n===========================\n") |
85 | ! |
cat("\nUser specified control parameters:") |
86 | ! |
cat("\n---------------------------\n") |
87 | ! |
for(i in x$usr.specified){ |
88 | ! |
cat(i, "=", x[[i]], fill = 80) |
89 |
} |
|
90 | ! |
cat("\nDefault control parameters:") |
91 | ! |
cat("\n---------------------------\n") |
92 | ! |
pToPrint <- names(x)[which(!grepl(".FUN", names(x)))] |
93 | ! |
pToPrint<- setdiff(pToPrint, c(x$usr.specified, "usr.specified")) |
94 | ! |
for (i in pToPrint) { |
95 | ! |
cat(i, "=", x[[i]], fill = 80) |
96 |
} |
|
97 | ! |
cat("Base Modules:", x$bi.mods, fill = 80) |
98 | ! |
if (length(x$user.mods) > 0) { |
99 | ! |
cat("Extension Modules:", x$user.mods, fill = 80) |
100 |
} |
|
101 |
|
|
102 | ! |
invisible() |
103 |
} |
|
104 | ||
105 |
#' @export |
|
106 |
print.seiqhrf <- function(x, ...) { |
|
107 |
|
|
108 | ! |
cat("SEIQHRF Model Simulation") |
109 | ! |
cat("\n=======================") |
110 | ! |
cat("\nModel class:", class(x)) |
111 |
|
|
112 | ! |
cat("\n\nSimulation Summary") |
113 | ! |
cat("\n-----------------------") |
114 | ! |
cat("\nModel type:", x$control$type) |
115 | ! |
cat("\nNumber of simulations:", x$control$nsims) |
116 | ! |
cat("\nNumber of time steps:", x$control$nsteps) |
117 |
|
|
118 | ! |
cat("\n \nModel Output (variable names)") |
119 | ! |
cat("\n-----------------------\n") |
120 | ! |
cat(names(x$epi), fill = 60) |
121 |
|
|
122 | ! |
invisible() |
123 |
} |
|
124 | ||
125 |
#' @export |
|
126 |
print.summary.seiqhrf <- function(x, ...){ |
|
127 |
|
|
128 | ! |
cat("Peaks of SEIQHRF Model Simulation: ") |
129 | ! |
cat("\n===============================\n") |
130 |
|
|
131 | ! |
No_comps <- length(x) |
132 | ! |
ret <- matrix(NA, No_comps, 2) |
133 | ! |
for(i in 1:nrow(ret)){ |
134 | ! |
tmp_mean <- x[[i]]$mean |
135 | ! |
tmp_mean[is.na(tmp_mean)] <- 0 |
136 | ! |
ret[i, 1] <- max(tmp_mean) |
137 | ! |
ret[i, 2] <- which(x[[i]]$mean == ret[i, 1])[1] |
138 |
} |
|
139 | ! |
colnames(ret) <- c("Max", "Time") |
140 | ! |
rownames(ret) <- names(x) |
141 | ||
142 | ! |
ret[, 1] <- format(ret[, 1], nsmall = 2, digits = 2) |
143 | ! |
print(ret, quote = FALSE, right = TRUE) |
144 |
} |
1 |
#' @title Cross Checking of Inputs for Stochastic Individual Contact Models |
|
2 |
#' |
|
3 |
#' @description This function checks that the three parameter lists from |
|
4 |
#' \code{\link{param_seiqhrf}}, \code{\link{init.icm}}, and |
|
5 |
#' \code{\link{control_seiqhrf}} are consistent. |
|
6 |
#' |
|
7 |
#' @param param An \code{EpiModel} object of class \code{param.seiqhrf}. |
|
8 |
#' @param init An \code{EpiModel} object of class \code{init.icm}. |
|
9 |
#' @param control A list returned by \code{\link{control_seiqhrf}}. |
|
10 |
#' |
|
11 |
#' @return |
|
12 |
#' This function returns no objects. |
|
13 |
#' |
|
14 |
#' @keywords internal |
|
15 |
#' |
|
16 |
crosscheck_seiqhrf <- function(param, init, control) { |
|
17 |
|
|
18 |
### Main class check |
|
19 | ! |
if (!inherits(param, "param.seiqhrf")) stop("param must an object of class param.seiqhrf", call. = FALSE) |
20 | ! |
if (!inherits(control, "control.seiqhrf")) stop("control must an object of class control.seiqhrf", call. = FALSE) |
21 |
|
|
22 |
|
|
23 |
|
|
24 | 2x |
if (control$skip.check == FALSE) { |
25 |
|
|
26 |
## Check that rec.rate is supplied for SIR models |
|
27 | 2x |
if (control$type %in% c("SIR", "SIS") && is.null(param$rec.rate)) { |
28 | ! |
stop("Specify rec.rate in param.seiqhrf", call. = FALSE) |
29 |
} |
|
30 |
|
|
31 |
|
|
32 |
## Check that paramets and init are supplied for SIR models |
|
33 | 2x |
if (control$type == "SIR" && is.null(init$r.num)) { |
34 | ! |
stop("Specify r.num in init.icm", call. = FALSE) |
35 |
} |
|
36 |
|
|
37 |
|
|
38 |
## Deprecated parameters |
|
39 | 2x |
bim <- grep(".FUN", names(control), value = TRUE) |
40 | 2x |
um <- which(grepl(".FUN", names(control)) & !(names(control) %in% bim)) |
41 | 2x |
if (length(um) == 0 && !is.null(control$type)) { |
42 | 2x |
if (!is.null(param$trans.rate)) { |
43 | ! |
stop("The trans.rate parameter is deprecated. Use the inf.prob ", |
44 | ! |
"parameter instead.", call. = FALSE) |
45 |
} |
|
46 |
} |
|
47 |
|
|
48 |
## Check param (modified from infection.FUN and arrivals.FUN) ------------------------------------- |
|
49 | 2x |
check_list <- c("act.rate.i", "inf.prob.i", "a.rate", "a.prop.e", "a.prop.i", "a.prop.q") |
50 | 2x |
if(control$type %in% c("SEIQHR", "SEIQHRF")) check_list<- c(check_list, "quar.rate", "disch.rate", "act.rate.e", "inf.prob.e", "act.rate.q", "inf.prob.q") |
51 | 2x |
check_idx <- which(names(param) %in% check_list) |
52 |
|
|
53 | 2x |
for(i in check_idx){ |
54 | 24x |
param_length = length(param[[i]]) |
55 | ! |
if (param_length != 1 && param_length != control$nsteps) stop(paste0("Length of", names(param)[i], "must be 1 or the value of nsteps")) |
56 |
} |
|
57 |
} |
|
58 |
|
|
59 |
## Check flare parameters |
|
60 | 2x |
if(!is.null(param$flare.inf.point)){ |
61 | ! |
if(length(param$flare.inf.point) != length(param$flare.inf.num)) stop("flare.inf.point and flare.inf.num must have the same length!") |
62 | ! |
if(max(param$flare.inf.point) > control$nsteps) stop("flare.inf.point can not exceed the total time points: nsteps!") |
63 | ! |
if(min(param$flare.inf.point) <= 2) stop("flare.inf.point should all be larger than 2, please change the initial i.num if flare time point needs to be smaller than 2") |
64 |
|
|
65 |
} |
|
66 |
## In-place assignment to update param and control |
|
67 | 2x |
on.exit(assign("param", param, pos = parent.frame())) |
68 | 2x |
on.exit(assign("control", control, pos = parent.frame()), add = TRUE) |
69 |
} |
1 |
#' Plot simulation result |
|
2 |
#' |
|
3 |
#' @param x An seiqhrf object returned from function \code{\link{seiqhrf}}. |
|
4 |
#' @param method If "times", plot Duration frequency distributions. |
|
5 |
#' If "weekly_local", plot local weekly estimates from simulation. |
|
6 |
#' If NULL, plot sirplus plots. |
|
7 |
#' @param return_df In effect only when method == "weekly", if TRUE returns |
|
8 |
#' also the dataframe used for plotting as well as the ggplot object. |
|
9 |
#' @param comp_remove Compartments to remove. Suggest c(s.num, r.num) |
|
10 |
#' @param time_limit Number of steps (days) to plot. |
|
11 |
#' @param ci T/F to include 95\% confidence intervals in sirplus plot. |
|
12 |
#' @param sep_compartments T/F use faceting to show each compartment in a |
|
13 |
#' separate plot, only works if plotting a single simulation. |
|
14 |
#' @param trans Y-axis transformation (e.g. log2, log10). |
|
15 |
#' @param known Dataframe with known compartment numbers to plot alongside |
|
16 |
#' projections |
|
17 |
#' @param start_date Date for day 0. |
|
18 |
#' @param show_start_date First date to show in plots. Use ymd format. If FALSE, |
|
19 |
#' shows from step 1. |
|
20 |
#' @param x_axis Title for x-axis. |
|
21 |
#' @param plot_title Title for whole plot. |
|
22 |
#' @param market.share between 0 and 1, percentage of local hospital beds in |
|
23 |
#' the simulated unit (e.g. state) |
|
24 |
#' @param icu_percent between 0 and 1, percentage of patients that should go to |
|
25 |
#' ICU among the ones that need hospitalization |
|
26 |
#' @param sim_population Size of population simulated. Only needed if providing |
|
27 |
#' `total_population`. |
|
28 |
#' @param total_population True population size, needed only if simulation size |
|
29 |
#' is smaller than the true population size due to computational cost |
|
30 |
#' etc. |
|
31 |
#' @param ... Additional parameters |
|
32 |
#' |
|
33 |
#' @importFrom tidyr pivot_longer |
|
34 |
#' @importFrom dplyr mutate |
|
35 |
#' @importFrom dplyr "%>%" |
|
36 |
#' @importFrom dplyr bind_rows |
|
37 |
#' @importFrom dplyr filter |
|
38 |
#' @importFrom lubridate ymd |
|
39 |
#' @import ggplot2 |
|
40 |
#' @export |
|
41 |
plot.seiqhrf <- function(x, method = NULL, |
|
42 |
comp_remove = "none", |
|
43 |
time_limit = 90, |
|
44 |
ci = TRUE, |
|
45 |
sep_compartments = FALSE, |
|
46 |
trans = FALSE, |
|
47 |
known = NULL, |
|
48 |
start_date = ymd("2020-03-21"), |
|
49 |
show_start_date = FALSE, |
|
50 |
x_axis = 'Date (MM-DD)', |
|
51 |
plot_title = '', |
|
52 |
return_df = TRUE, |
|
53 |
market.share = .04, |
|
54 |
icu_percent = .1, |
|
55 |
sim_population = 1000, |
|
56 |
total_population = NULL, ...) { |
|
57 |
|
|
58 | ! |
if(is.null(method)){ |
59 | ! |
plot_sirplus(x, comp_remove = comp_remove, |
60 | ! |
time_limit = time_limit, |
61 | ! |
ci = ci, |
62 | ! |
sep_compartments = sep_compartments, |
63 | ! |
trans = trans, |
64 | ! |
known = known, |
65 | ! |
start_date = start_date, |
66 | ! |
x_axis = x_axis, |
67 | ! |
plot_title = plot_title, |
68 | ! |
sim_population = sim_population, |
69 | ! |
total_population = total_population) |
70 |
|
|
71 | ! |
}else if(method == "times"){ |
72 |
|
|
73 | ! |
if(!inherits(x, "seiqhrf")) stop( |
74 | ! |
"If method == times, x needs to be an seiqhrf object") |
75 | ! |
plot_times(x) |
76 |
|
|
77 | ! |
}else if(method == "weekly_local"){ |
78 |
|
|
79 | ! |
ret <- get_weekly_local(x, market.share = market.share, |
80 | ! |
icu_percent = icu_percent, |
81 | ! |
start_date = start_date, |
82 | ! |
show_start_date = show_start_date, |
83 | ! |
time_limit = time_limit, |
84 | ! |
sim_population = sim_population, |
85 | ! |
total_population = total_population) |
86 | ! |
if(return_df){ |
87 | ! |
return(ret) |
88 |
}else{ |
|
89 | ! |
return(ret$plot) |
90 |
} |
|
91 |
} |
|
92 |
} |
|
93 | ||
94 |
#' Plot simulation result given list |
|
95 |
#' |
|
96 |
#' @param x A list of seiqhrf objects returned from \code{\link{seiqhrf}}. |
|
97 |
#' @param comp_remove Compartments to remove. Suggest c(s.num, r.num) |
|
98 |
#' @param time_limit Number of steps (days) to plot. |
|
99 |
#' @param ci T/F to include 95\% confidence intervals in sirplus plot. |
|
100 |
#' @param sep_compartments T/F use faceting to show each compartment in a |
|
101 |
#' separate plot, only works if plotting a single simulation. |
|
102 |
#' @param trans Y-axis transformation (e.g. log2, log10). |
|
103 |
#' @param known Dataframe with known compartment numbers to plot alongside |
|
104 |
#' projections |
|
105 |
#' @param start_date Date for day 0. |
|
106 |
#' @param x_axis Title for x-axis. |
|
107 |
#' @param plot_title Title for whole plot. |
|
108 |
#' total_population |
|
109 |
#' @param sim_population Size of population simulated. Only needed if providing |
|
110 |
#' `total_population`. |
|
111 |
#' @param total_population True population size, needed only if simulation size |
|
112 |
#' is smaller than the true population size due to computational cost |
|
113 |
#' etc. |
|
114 |
#' @param ... Additional parameters |
|
115 |
#' |
|
116 |
#' @importFrom tidyr pivot_longer |
|
117 |
#' @importFrom dplyr mutate |
|
118 |
#' @importFrom dplyr "%>%" |
|
119 |
#' @importFrom dplyr bind_rows |
|
120 |
#' @importFrom dplyr filter |
|
121 |
#' @importFrom lubridate ymd |
|
122 |
#' @import ggplot2 |
|
123 |
#' @export |
|
124 |
plot.list <- function(x, comp_remove = "none", |
|
125 |
time_limit = 90, |
|
126 |
ci = TRUE, |
|
127 |
sep_compartments = FALSE, |
|
128 |
trans = FALSE, |
|
129 |
known = NULL, |
|
130 |
start_date = ymd("2020-03-21"), |
|
131 |
x_axis = 'Date (MM-DD)', |
|
132 |
plot_title = '', |
|
133 |
sim_population = 1000, |
|
134 |
total_population = NULL, ...){ |
|
135 |
|
|
136 | 1x |
plot_sirplus(x, comp_remove = comp_remove, |
137 | 1x |
time_limit = time_limit, |
138 | 1x |
ci = ci, |
139 | 1x |
sep_compartments = sep_compartments, |
140 | 1x |
trans = trans, |
141 | 1x |
known = known, |
142 | 1x |
start_date = start_date, |
143 | 1x |
x_axis = x_axis, |
144 | 1x |
plot_title = plot_title, |
145 | 1x |
sim_population = sim_population, |
146 | 1x |
total_population = total_population) |
147 |
} |
|
148 | ||
149 | ||
150 | ||
151 | ||
152 |
#' Wrapper for primary sirplus plotting function |
|
153 |
#' |
|
154 |
#' Flexible function to generate sirplus plots (i.e. compartment counts over |
|
155 |
#' time). This function allows for plotting multiple experiments, viewing the |
|
156 |
#' plots of different scales (e.g. log2), plotting compartments separately, |
|
157 |
#' adding 95\% CIs, and plotting known data along side the simulations. |
|
158 |
#' |
|
159 |
#' @param x A seiqhrf object (or list of multiple seiqhrf objects) returned |
|
160 |
#' from \code{\link{seiqhrf}}. |
|
161 |
#' @param comp_remove Compartments to remove. Suggest c(s.num, r.num) |
|
162 |
#' @param time_limit Number of steps (days) to plot. |
|
163 |
#' @param ci T/F to include 95\% confidence intervals in sirplus plot. |
|
164 |
#' @param sep_compartments T/F use faceting to show each compartment in a |
|
165 |
#' separate plot, only works if plotting a single simulation. |
|
166 |
#' @param trans Y-axis transformation (e.g. log2, log10). |
|
167 |
#' @param known Dataframe with known compartment numbers to plot alongside |
|
168 |
#' projections |
|
169 |
#' @param start_date Date for day 0. |
|
170 |
#' @param x_axis Title for x-axis. |
|
171 |
#' @param plot_title Title for whole plot. |
|
172 |
#' @param sim_population Size of population simulated. Only needed if providing |
|
173 |
#' `total_population`. |
|
174 |
#' @param total_population True population size, needed only if simulation size |
|
175 |
#' is smaller than the true population size due to computational cost |
|
176 |
#' @param ... Additional parameters |
|
177 |
#' |
|
178 |
#' @return ggplot2 object |
|
179 |
#' |
|
180 |
#' @importFrom tidyr pivot_longer |
|
181 |
#' @importFrom dplyr mutate |
|
182 |
#' @importFrom dplyr "%>%" |
|
183 |
#' @importFrom dplyr bind_rows |
|
184 |
#' @importFrom dplyr filter |
|
185 |
#' @import ggplot2 |
|
186 |
#' @export |
|
187 |
plot_sirplus <- function(x, comp_remove, |
|
188 |
time_limit, |
|
189 |
ci, |
|
190 |
sep_compartments, |
|
191 |
trans, |
|
192 |
known, |
|
193 |
start_date, |
|
194 |
x_axis, |
|
195 |
plot_title, |
|
196 |
sim_population, |
|
197 |
total_population,...){ |
|
198 |
|
|
199 |
# Convert from seiqhrf object to dataframe |
|
200 | 1x |
plot_df <- format_sims(x, time_limit = time_limit, start_date = start_date) |
201 | ! |
reo_exp <- function(x) {factor(x, levels = unique(plot_df$experiment))} |
202 |
|
|
203 |
# Get Confidence Intervals |
|
204 | 1x |
if(ci){ |
205 | 1x |
plot_df <- get_ci(x, plot_df) |
206 |
} |
|
207 |
|
|
208 |
# Scale up to full population size if needed |
|
209 | 1x |
if(!is.null(total_population)){ |
210 | ! |
if(total_population < sim_population) |
211 | ! |
stop("total population should be larger than simulated size") |
212 |
|
|
213 | ! |
scale_factor <- total_population/sim_population |
214 | ! |
plot_df <- plot_df %>% mutate(count = count * scale_factor) |
215 |
|
|
216 | ! |
if(ci){ |
217 | ! |
plot_df <- plot_df %>% mutate(CI.1 = count * scale_factor, |
218 | ! |
CI.2 = CI.2 * scale_factor, |
219 | ! |
qntCI.1 = qntCI.1 * scale_factor, |
220 | ! |
qntCI.2 = qntCI.2 * scale_factor) |
221 |
} |
|
222 |
} |
|
223 |
|
|
224 |
# Add known compartment counts |
|
225 | 1x |
if(is.data.frame(known)){ |
226 | ! |
plot_df <- add_known(plot_df, known = known, start_date = start_date) |
227 |
} |
|
228 |
|
|
229 |
# Define compartment names and colours |
|
230 | 1x |
comps <- c("s.num", "e.num", "i.num", "q.num", "h.num", "r.num", "f.num") |
231 | 1x |
compcols <- c(s.num = "#4477AA", e.num = "#66CCEE", i.num = "#CCBB44", |
232 | 1x |
q.num = "#AA3377", h.num = "#EE6677", r.num = "#228833", |
233 | 1x |
f.num = "#BBBBBB") |
234 | 1x |
complabels <- c(s.num = "S: Susceptible", e.num = "E: Asymptomatic", |
235 | 1x |
i.num = "I: Infected", q.num = "Q: Self-isolated", |
236 | 1x |
h.num = "H: Hospitalized", r.num = "R: Recovered", |
237 | 1x |
f.num = "F: Case Fatalities") |
238 |
|
|
239 |
# Filter compartments |
|
240 | 1x |
comp_plot <- setdiff(comps, comp_remove) |
241 | 1x |
plot_df <- plot_df %>% filter(compartment %in% c(comp_plot)) %>% |
242 | 1x |
filter(time <= time_limit) |
243 |
|
|
244 |
# Plot with options |
|
245 | 1x |
p <- ggplot(plot_df, aes(x = Date, y = count, colour = compartment, |
246 | 1x |
linetype = sim)) + |
247 | 1x |
geom_line(size = 1.2, alpha = 0.8) + |
248 | 1x |
scale_x_date(date_breaks = "1 week", date_labels = "%m-%d") + |
249 | 1x |
scale_colour_manual(values = compcols, labels = complabels) + |
250 | 1x |
labs(title = plot_title, x = x_axis, y = "Prevalence") + |
251 | 1x |
theme_bw() + theme(axis.text.x = element_text(angle = 90)) |
252 |
|
|
253 | 1x |
if(length(unique(plot_df$experiment)) > 1){ |
254 | 1x |
p <- p + facet_grid(reo_exp(experiment) ~ ., scales = 'free') |
255 |
} |
|
256 |
|
|
257 | 1x |
if(sep_compartments){ |
258 | ! |
p <- p + facet_grid(compartment ~ ., scales = 'free') |
259 |
} |
|
260 |
|
|
261 | 1x |
if(trans){ |
262 | ! |
p <- p + scale_y_continuous(trans = trans) |
263 |
} |
|
264 |
|
|
265 | 1x |
if(ci){ |
266 | 1x |
p <- p + geom_ribbon(aes(ymin=qntCI.1, ymax=qntCI.2, x=Date, |
267 | 1x |
fill = compartment, colour = NULL), |
268 | 1x |
alpha = 0.4) + |
269 | 1x |
scale_color_manual(values = compcols, labels = complabels) + |
270 | 1x |
scale_fill_manual(values = compcols, guide = FALSE) |
271 |
} |
|
272 |
|
|
273 | 1x |
p |
274 |
} |
|
275 | ||
276 |
#' Plot compartment duration distributions |
|
277 |
#' |
|
278 |
#' Function to plot Duration frequency distributions. If multiple simulations |
|
279 |
#' were performed (nsim >1), durations from sims are appended to each other. |
|
280 |
#' |
|
281 |
#' @param sim An seiqhrf object returned from function \code{\link{seiqhrf}}. |
|
282 |
#' |
|
283 |
#' @return ggplot2 object |
|
284 |
#' |
|
285 |
#' @import ggplot2 |
|
286 |
#' @importFrom tidyr pivot_longer |
|
287 |
#' @importFrom dplyr mutate |
|
288 |
#' @importFrom dplyr "%>%" |
|
289 |
#' @importFrom dplyr bind_rows |
|
290 |
#' @importFrom dplyr select |
|
291 |
#' @importFrom dplyr filter |
|
292 |
#' |
|
293 |
#' @export |
|
294 |
plot_times <- function(sim) { |
|
295 |
|
|
296 | ! |
for (s in 1:sim$control$nsims) { |
297 | ! |
if (s == 1) { |
298 | ! |
times <- sim$times[[paste("sim", s, sep = "")]] |
299 | ! |
times <- times %>% mutate(s = s) |
300 |
} else { |
|
301 | ! |
times <- times %>% bind_rows(sim$times[[paste("sim", s, sep = "")]] |
302 | ! |
%>% mutate(s = s)) |
303 |
} |
|
304 |
} |
|
305 |
|
|
306 | ! |
times <- times %>% mutate(infTime = ifelse(infTime < 0, -5, infTime), |
307 | ! |
expTime = ifelse(expTime < 0, -5, expTime)) %>% |
308 | ! |
mutate(incubation_period = infTime - expTime, |
309 | ! |
illness_duration = recovTime - expTime, |
310 | ! |
illness_duration_hosp = dischTime - expTime, |
311 | ! |
hosp_los = dischTime - hospTime, |
312 | ! |
quarantine_delay = quarTime - infTime, |
313 | ! |
survival_time = fatTime - infTime) %>% |
314 | ! |
select(s, incubation_period, quarantine_delay, illness_duration, |
315 | ! |
illness_duration_hosp, hosp_los, survival_time) %>% |
316 | ! |
pivot_longer(-s, names_to = "period_type", values_to = "duration") %>% |
317 | ! |
mutate(period_type = factor(period_type, |
318 | ! |
levels = c("incubation_period", |
319 | ! |
"quarantine_delay", |
320 | ! |
"illness_duration", |
321 | ! |
"illness_duration_hosp", |
322 | ! |
"hosp_los", |
323 | ! |
"survival_time"), |
324 | ! |
labels = c("Incubation\nperiod", |
325 | ! |
"Delay entering\nisolation", |
326 | ! |
"Illness\nduration", |
327 | ! |
"Illness\nduration (hosp)", |
328 | ! |
"Hospital stay\nduration", |
329 | ! |
"Survival time\nfor fatalities"), |
330 | ! |
ordered = TRUE)) |
331 |
|
|
332 | ! |
times %>% filter(duration <= 30) %>% ggplot(aes(x = duration)) + |
333 | ! |
geom_bar() + facet_grid(period_type ~ ., scales = "free_y") + |
334 | ! |
labs(title = "Compartment Duration Distributions") |
335 |
} |
|
336 | ||
337 | ||
338 |
#' Extract and plot information of local and weekly estimates from simulation |
|
339 |
#' |
|
340 |
#' @param sim An \code{seiqhrf} object returned by \link{simulate_seiqhrf}. |
|
341 |
#' @param market.share between 0 and 1, percentage of local hospital beds in |
|
342 |
#' the simulated unit (e.g. state) |
|
343 |
#' @param icu_percent between 0 and 1, percentage of patients that should go to |
|
344 |
#' ICU among the ones that need hospitalization |
|
345 |
#' @param start_date Epidemic start date. Default is 'na', if not provided will |
|
346 |
#' plot week numbers, if provided will plot the first day (Sunday) of the |
|
347 |
#' week. |
|
348 |
#' @param show_start_date First date to show in plots. Use ymd format. If FALSE, |
|
349 |
#' shows from step 1. |
|
350 |
#' @param time_limit Number of days to include. |
|
351 |
#' @param sim_population Size of population simulated. Only needed if providing |
|
352 |
#' `total_population`. |
|
353 |
#' @param total_population True population size, needed only if simulation size |
|
354 |
#' is smaller than the true population size due to computational cost |
|
355 |
#' etc. |
|
356 |
#' |
|
357 |
#' @return |
|
358 |
#' \itemize{ |
|
359 |
#' \item \code{plot:} A \code{ggplot} object, bar charts of count of patients |
|
360 |
#' requiring hospitalization and ICU respectively |
|
361 |
#' \item \code{result:} A dataframe |
|
362 |
#' \itemize{\item \code{week:} week number from input \code{sim}, |
|
363 |
#' \item \code{hosp:} the number of patients that require hospitalization locally, |
|
364 |
#' \item \code{icu:} the number of patients that require ICU locally. } |
|
365 |
# |
|
366 |
#' } |
|
367 |
#' |
|
368 |
#' @importFrom tidyr pivot_wider |
|
369 |
#' @importFrom dplyr group_by |
|
370 |
#' @importFrom dplyr summarise |
|
371 |
#' @importFrom utils head |
|
372 |
#' @export |
|
373 |
get_weekly_local <- function(sim, |
|
374 |
market.share, |
|
375 |
icu_percent, |
|
376 |
start_date, |
|
377 |
show_start_date, |
|
378 |
time_limit, |
|
379 |
sim_population, |
|
380 |
total_population){ |
|
381 | ||
382 |
# Get h.num and 95% quantile CIs |
|
383 | ! |
sim_mean <- as.data.frame(sim, out = "mean") |
384 | ! |
ci_info <- as.data.frame.list(summary.seiqhrf(sim)$h.num) |
385 | ! |
hosp <- data.frame('h.num' = sim_mean$h.num, 'ci5' = ci_info$qntCI.1, |
386 | ! |
'ci95' = ci_info$qntCI.2) |
387 | ! |
hosp[is.na(hosp)] <- 0 |
388 | ! |
hosp <- hosp[1: time_limit, ] |
389 | ! |
hosp$date <- start_date + as.numeric(row.names(hosp)) |
390 |
|
|
391 | ! |
if(class(show_start_date) == "Date"){ |
392 | ! |
hosp <- hosp %>% filter(date >= show_start_date) |
393 |
} |
|
394 |
|
|
395 |
# Scale for population size and hospital market share if needed |
|
396 | ! |
if(!is.null(total_population)){ |
397 | ! |
if(total_population < sim_population) |
398 | ! |
stop("total Population should be larger than simulated size") |
399 | ! |
cat("Scalling w.r.t total population") |
400 | ! |
date_tmp <- hosp$date |
401 | ! |
hosp$date <- NULL |
402 | ! |
print(utils::head(hosp)) |
403 | ! |
hosp <- hosp*total_population/sim_population |
404 | ! |
hosp$date <- date_tmp |
405 |
} |
|
406 |
|
|
407 | ! |
if(market.share < 0 || market.share > 1) stop("Market share has to be between |
408 | ! |
0 and 1") |
409 | ! |
if(icu_percent < 0 || icu_percent > 1) stop("ICU percentage has to be between |
410 | ! |
0 and 1") |
411 |
|
|
412 |
# Get weekly sums & calculate projected icu numbers |
|
413 | ! |
hosp.wk <- hosp %>% group_by(yr_wk = cut(date, "week", start.on.monday = TRUE)) %>% |
414 | ! |
summarise(h.num = sum(h.num), h.ci5 = sum(ci5), h.ci95=sum(ci95)) %>% |
415 | ! |
mutate(icu.num = h.num * icu_percent, |
416 | ! |
icu.ci5 = h.ci5 * icu_percent, |
417 | ! |
icu.ci95 = h.ci95 * icu_percent) %>% |
418 | ! |
mutate(h.num = h.num - icu.num, |
419 | ! |
h.ci5 = h.ci5 - icu.ci5, |
420 | ! |
h.ci95 = h.ci95 - icu.ci95) |
421 |
|
|
422 |
# Make long format for ggplot |
|
423 | ! |
hosp.wk2 <- hosp.wk %>% pivot_longer(-yr_wk, names_to = c('type', 'metric'), |
424 | ! |
values_to = 'val', |
425 | ! |
names_pattern = '(h|icu).(num|ci5|ci95)') %>% |
426 | ! |
pivot_wider(names_from = metric, values_from = val) |
427 | ||
428 | ! |
p <- ggplot(hosp.wk2, aes(x = yr_wk, y = num, color = type)) + |
429 | ! |
geom_point() + |
430 | ! |
labs(y="Weekly Cumulative Count", x = "Week Start (Monday: YYYY-MM-DD)") + |
431 | ! |
geom_errorbar(aes(ymin=ci5, ymax=ci95), size=0.5, width=.25) + |
432 | ! |
scale_color_discrete(name = "Type", labels = c("General", "ICU")) + |
433 | ! |
theme_bw() + theme(axis.text.x = element_text(angle = 90)) |
434 |
|
|
435 | ! |
return(list("plot" = p, "result" = hosp.wk)) |
436 |
|
|
437 |
} |
|
438 | ||
439 | ||
440 | ||
441 |
#' Format seiqhrf objects into dataframe for ggplot |
|
442 |
#' |
|
443 |
#' @param x An seiqhrf object returned from function \code{\link{seiqhrf}}. |
|
444 |
#' @param time_limit Number of steps (days) to plot. |
|
445 |
#' @param start_date Date for day 0. |
|
446 |
#' |
|
447 |
#' @return dataframe |
|
448 |
#' |
|
449 |
#' @export |
|
450 |
format_sims <- function(x, time_limit, start_date){ |
|
451 |
|
|
452 |
# Merge models to plot together |
|
453 | 1x |
if(class(x) == "seiqhrf"){ |
454 | ! |
sim_id <- "seiqhrf model" |
455 | ! |
plot_df <- as.data.frame(x, out = "mean") |
456 | ! |
plot_df <- plot_df %>% mutate(experiment = sim_id) |
457 |
}else{ |
|
458 |
|
|
459 | 1x |
sim_id <- names(x) |
460 | ! |
if(is.null(sim_id)) stop("Please assign a name to each element in sims") |
461 |
|
|
462 | 1x |
plot_df <- as.data.frame(x[[1]], out = "mean") |
463 | 1x |
plot_df <- plot_df %>% mutate(experiment = sim_id[1]) |
464 | 1x |
if(length(sim_id) > 1){ |
465 | 1x |
for (i in (2:length(sim_id))) { |
466 | 1x |
tmp_df <- as.data.frame(x[[i]], out = "mean") |
467 | 1x |
tmp_df <- tmp_df %>% mutate(experiment = sim_id[i]) |
468 | 1x |
plot_df <- plot_df %>% bind_rows(plot_df, tmp_df) |
469 |
} |
|
470 |
} |
|
471 |
} |
|
472 |
|
|
473 | 1x |
plot_df <- plot_df %>% filter(time <= time_limit) %>% |
474 | 1x |
pivot_longer(-c(time, experiment), |
475 | 1x |
names_to = "compartment", |
476 | 1x |
values_to = "count", |
477 | 1x |
values_ptypes = list(compartment = 'character', |
478 | 1x |
count = numeric())) %>% |
479 | 1x |
mutate(sim = 'sim', |
480 | 1x |
Date = start_date + time) |
481 | ||
482 | 1x |
return(plot_df) |
483 |
} |
|
484 | ||
485 |
#' Internal |
|
486 |
#' @param obj An seiqhrf object |
|
487 |
#' @param exp_name Name of experiment |
|
488 | ||
489 |
ci_info_update <- function(obj, exp_name = "seiqhrf model"){ |
|
490 |
|
|
491 | 2x |
ci_tmp <- as.data.frame.list(summary(obj)) |
492 | 2x |
ci_tmp <- ci_tmp %>% |
493 | 2x |
mutate(time = as.numeric(row.names(ci_tmp))) %>% |
494 | 2x |
pivot_longer(cols = -time, names_to = 'compartment', |
495 | 2x |
values_to = 'mean', |
496 | 2x |
values_ptypes = list(compartment = 'character', |
497 | 2x |
mean = numeric())) %>% |
498 | 2x |
tidyr::separate(compartment, |
499 | 2x |
into = c('compartment', 'metric'), |
500 | 2x |
sep='num.') %>% |
501 | 2x |
mutate(compartment = paste0(compartment, 'num'), |
502 | 2x |
experiment = exp_name) %>% |
503 | 2x |
pivot_wider(names_from = metric, values_from = mean) |
504 |
|
|
505 | 2x |
ci_tmp |
506 |
} |
|
507 | ||
508 |
#' Get 95\% confidence intervals |
|
509 |
#' |
|
510 |
#' @param x An seiqhrf object returned from function \code{\link{seiqhrf}} or a list of seiqhrf objects. |
|
511 |
#' @param plot_df Dataframe with known compartment numbers to plot alongside |
|
512 |
#' projections |
|
513 |
#' |
|
514 |
#' @return dataframe with CIs and sd added |
|
515 |
#' @importFrom tidyr separate |
|
516 |
#' @importFrom dplyr full_join |
|
517 |
#' |
|
518 |
#' @export |
|
519 |
#' |
|
520 |
get_ci <- function(x, plot_df){ |
|
521 |
|
|
522 |
# Get sim variance metrics for single seiqhrf object |
|
523 | 1x |
if(class(x) == "seiqhrf"){ |
524 | ! |
ci_info <- ci_info_update(x) |
525 |
}else{ |
|
526 |
|
|
527 | 1x |
sim_id <- names(x) |
528 | ! |
if(class(x) != "list") stop("The class of x should either be list or seiqhrf") |
529 | 1x |
ci_info <- ci_info_update(x[[1]], sim_id[1]) |
530 |
|
|
531 | 1x |
if(length(sim_id) > 1){ |
532 | 1x |
for (i in (2:length(sim_id))) { |
533 | 1x |
ci_tmp <- ci_info_update(x[[i]], sim_id[i]) |
534 | 1x |
ci_info <- ci_info %>% bind_rows(ci_info, ci_tmp) |
535 |
} |
|
536 |
} |
|
537 |
} |
|
538 |
|
|
539 | 1x |
ci_info <- ci_info %>% mutate(sim = 'sim') %>% |
540 | 1x |
mutate(sim = 'sim') |
541 | 1x |
ci_info[is.na(ci_info)] <- 0 |
542 |
|
|
543 | 1x |
plot_df <- plot_df %>% |
544 | 1x |
dplyr::full_join(ci_info, by = c('time', 'compartment', 'experiment', 'sim')) |
545 |
|
|
546 | 1x |
return(plot_df) |
547 |
} |
|
548 | ||
549 |
#' Add known counts to sims dataframe for ggplot |
|
550 |
#' |
|
551 |
#' @param plot_df An seiqhrf object returned from function \code{\link{seiqhrf}}. |
|
552 |
#' @param known Dataframe with known compartment numbers to plot alongside |
|
553 |
#' projections |
|
554 |
#' @param start_date Date for day 0. |
|
555 |
#' |
|
556 |
#' @return dataframe with known data added |
|
557 |
#' |
|
558 |
#' @export |
|
559 |
#' |
|
560 |
add_known <- function(plot_df, known, start_date){ |
|
561 |
|
|
562 |
# Add Date to known data |
|
563 | ! |
missing_cols <- setdiff(names(plot_df), names(known)) |
564 |
|
|
565 | ! |
if("Date" %in% missing_cols){ |
566 | ! |
known$Date = start_date + known$time |
567 |
} |
|
568 |
|
|
569 | ! |
known <- known %>% pivot_longer(cols = -c(time, Date), |
570 | ! |
names_to = 'compartment', |
571 | ! |
values_to = 'count', |
572 | ! |
values_ptypes = list(compartment = 'character', |
573 | ! |
count = numeric())) |
574 | ! |
exps <- unique(plot_df$experiment) |
575 | ! |
for(i in exps){ |
576 | ! |
known <- known %>% mutate(experiment = i, sim = 'known') |
577 | ! |
missing_cols <- setdiff(names(plot_df), names(known)) |
578 | ||
579 | ! |
if(length(missing_cols) > 0){ |
580 | ! |
add <- rep(0, length(missing_cols)) |
581 | ! |
names(add) <- missing_cols |
582 | ! |
known <- known %>% mutate(!!! add) |
583 |
} |
|
584 | ! |
plot_df <- rbind(plot_df, known) |
585 |
} |
|
586 |
|
|
587 | ! |
return(plot_df) |
588 |
} |
|
589 |
|
|
590 |
1 |
#' Plot models |
|
2 |
#' |
|
3 |
#' Function to plot individuals models or multiple models for comparison. |
|
4 |
#' |
|
5 |
#' @param sims Single or list of sims to plot |
|
6 |
#' @param sim_id String or list of strings to name each facet |
|
7 |
#' @param comp_remove Compartments to remove. Suggest c(s.num, r.num) |
|
8 |
#' @param time_lim Number of steps (days) to plot. |
|
9 |
#' @param trans Y-axis transformation (e.g. log2, log10). Default = none. |
|
10 |
#' @param known Dataframe with known compartment numbers to plot alongside |
|
11 |
#' projections |
|
12 |
#' @param start_date Date for day 0. Default: ymd("2020-03-21"), |
|
13 |
#' @param x_axis Title for x-axis. Default: 'Days since beginning of epidemic' |
|
14 |
#' @param plot_title Title for whole plot. Default: 'SEIQHRF plot' |
|
15 |
#' |
|
16 |
#' @return ggplot2 object |
|
17 |
#' |
|
18 |
#' @importFrom tidyr pivot_longer |
|
19 | ||
20 |
plot_models <- function(sims, |
|
21 |
sim_id = "baseline", |
|
22 |
comp_remove = "none", |
|
23 |
time_lim = 100, |
|
24 |
trans = 'na', |
|
25 |
known = NULL, |
|
26 |
start_date = ymd("2020-03-21"), |
|
27 |
x_axis = 'Days since beginning of epidemic', |
|
28 |
plot_title = 'SEIQHRF'){ |
|
29 |
|
|
30 |
# Define a standard set of colours to represent compartments |
|
31 | ! |
comps <- c("s.num", "e.num", "i.num", "q.num", "h.num", "r.num", "f.num") |
32 | ! |
compcols <- c(s.num = "#4477AA", e.num = "#66CCEE", i.num = "#CCBB44", |
33 | ! |
q.num = "#AA3377", h.num = "#EE6677", r.num = "#228833", |
34 | ! |
f.num = "#BBBBBB") |
35 | ! |
complabels <- c(s.num = "S: Susceptible", e.num = "E: Asymptomatic", |
36 | ! |
i.num = "I: Infected", q.num = "Q: Self-isolated", |
37 | ! |
h.num = "H: Hospitalized", r.num = "R: Recovered", |
38 | ! |
f.num = "F: Case Fatalities") |
39 |
|
|
40 |
# Merge models to plot together |
|
41 | ! |
for (i in (1: length(sim_id))) { |
42 | ! |
if (i == 1) { |
43 | ! |
plot_df <- sims[i*2]$df |
44 | ! |
plot_df <- plot_df %>% mutate(experiment = sim_id[i]) |
45 |
} else{ |
|
46 | ! |
tmp_df <- sims[i*2]$df |
47 | ! |
tmp_df <- tmp_df %>% mutate(experiment = sim_id[i]) |
48 | ! |
plot_df <- plot_df %>% bind_rows(plot_df, tmp_df) |
49 |
} |
|
50 |
} |
|
51 |
|
|
52 | ! |
plot_df <- plot_df %>% filter(time <= time_lim) %>% |
53 | ! |
pivot_longer(-c(time, experiment), |
54 | ! |
names_to = "compartment", |
55 | ! |
values_to = "count") |
56 | ! |
plot_df$sim <- 'sim' |
57 | ! |
plot_df$Date <- start_date + plot_df$time |
58 | ! |
plot_df$time <- NULL |
59 |
|
|
60 |
# Add known data |
|
61 | ! |
if(is.data.frame(known)){ |
62 | ! |
exps <- unique(plot_df$experiment) |
63 | ! |
for(i in exps){ |
64 | ! |
known$experiment <- i |
65 | ! |
plot_df <- rbind(plot_df, known) |
66 |
} |
|
67 |
} |
|
68 |
|
|
69 |
# Filter compartments |
|
70 | ! |
comp_plot <- setdiff(comps, comp_remove) |
71 | ! |
plot_df <- plot_df %>% filter(compartment %in% c(comp_plot)) |
72 | ||
73 | ! |
reo_exp <- function(x) { factor(x, levels = sim_id)} |
74 |
|
|
75 |
|
|
76 |
# Plot single model |
|
77 | ! |
if(trans == "na"){ |
78 | ! |
if(length(sim_id) == 1){ |
79 | ! |
plot_df %>% ggplot(aes(x = Date, y = count+1, colour = compartment, |
80 | ! |
linetype = sim)) + |
81 | ! |
geom_line(size = 1.5, alpha = 0.8) + |
82 | ! |
scale_x_date(date_breaks = "1 week", date_labels = "%m-%d") + |
83 | ! |
scale_colour_manual(values = compcols, labels = complabels) + |
84 | ! |
labs(title = plot_title, x = x_axis, y = "Prevalence") + |
85 | ! |
theme_bw() + |
86 | ! |
theme(axis.text.x = element_text(angle = 90)) |
87 |
} else ( |
|
88 | ! |
plot_df %>% ggplot(aes(x = Date, y = count, colour = compartment, |
89 | ! |
linetype = sim)) + |
90 | ! |
facet_grid(reo_exp(experiment) ~ ., scales = 'free') + |
91 | ! |
scale_x_date(date_breaks = "1 week", date_labels = "%m-%d") + |
92 | ! |
geom_line(size = 1.5, alpha = 0.8) + |
93 | ! |
scale_colour_manual(values = compcols, labels = complabels) + |
94 | ! |
labs(title = plot_title, x = x_axis, y = "Prevalence") + |
95 | ! |
theme_bw() + |
96 | ! |
theme(axis.text.x = element_text(angle = 90)) |
97 |
) |
|
98 |
}else{ |
|
99 | ! |
if(length(sim_id) == 1){ |
100 | ! |
plot_df %>% ggplot(aes(x = Date, y = count+1, colour = compartment, |
101 | ! |
linetype = sim)) + |
102 | ! |
scale_y_continuous(trans = trans) + |
103 | ! |
scale_x_date(date_breaks = "1 week", date_labels = "%m-%d") + |
104 | ! |
geom_line(size = 1.5, alpha = 0.8) + |
105 | ! |
scale_colour_manual(values = compcols, labels = complabels) + |
106 | ! |
labs(title = plot_title, x = x_axis, y = "Prevalence") + |
107 | ! |
theme_bw() + |
108 | ! |
theme(axis.text.x = element_text(angle = 90)) |
109 |
} else ( |
|
110 | ! |
plot_df %>% ggplot(aes(x = Date, y = count+1, colour = compartment, |
111 | ! |
linetype = sim)) + |
112 | ! |
facet_grid(reo_exp(experiment) ~ ., scales = 'free') + |
113 | ! |
scale_y_continuous(trans = trans) + |
114 | ! |
scale_x_date(date_breaks = "1 week", date_labels = "%m-%d") + |
115 | ! |
geom_line(size = 1.5, alpha = 0.8) + |
116 | ! |
scale_colour_manual(values = compcols, labels = complabels) + |
117 | ! |
labs(title = plot_title, x = x_axis, y = "Prevalence") + |
118 | ! |
theme_bw() + |
119 | ! |
theme(axis.text.x = element_text(angle = 90)) |
120 |
) |
|
121 |
} |
|
122 |
} |
1 |
#' Departures function |
|
2 |
#' |
|
3 |
#' Function to handle background demographics for the SEIQHRF model. |
|
4 |
#' Specifically departures (deaths not due to the virus, and emigration). |
|
5 |
#' |
|
6 |
#' @param dat Merged input parameters. |
|
7 |
#' @param at Step number |
|
8 |
#' @param seed random seed for checking consistency |
|
9 |
#' |
|
10 |
#' @return Updated dat |
|
11 |
#' @export |
|
12 |
departures.FUN <- function(dat, at, seed = NULL) { |
|
13 |
|
|
14 | 10x |
if(!is.null(seed)) set.seed(seed) |
15 |
|
|
16 |
# Conditions -------------------------------------------------------------- |
|
17 | ! |
if (!dat$param$vital) return(dat) |
18 |
|
|
19 |
# Variables ----------------------------------------------------------------- |
|
20 | 11x |
rate <- dat$param$ds.rate |
21 | 11x |
rand <- dat$control$d.rand |
22 | 11x |
status <- dat$attr$status |
23 | 11x |
active <- dat$attr$active |
24 |
|
|
25 |
# Susceptible departures ------------------------------------------------------ |
|
26 | 11x |
res <- update_active(rate, rand, active, status, label = "s") |
27 | 11x |
nDepartures <- res[[1]] |
28 | 1x |
if(!is.null(res[[2]])) active <- dat$attr$active[res[[2]]] <- 0 |
29 |
|
|
30 | ! |
if (at == 2) dat$epi$ds.flow <- c(0, nDepartures) else dat$epi$ds.flow[at] <- nDepartures |
31 |
|
|
32 |
# Exposed Departures --------------------------------------------------------- |
|
33 | 11x |
res <- update_active(rate, rand, active, status, label = "e") |
34 | 11x |
nDepartures <- res[[1]] |
35 | ! |
if(!is.null(res[[2]])) active <- dat$attr$active[res[[2]]] <- 0 |
36 |
|
|
37 | ! |
if (at == 2) dat$epi$de.flow <- c(0, nDepartures) else dat$epi$de.flow[at] <- nDepartures |
38 |
|
|
39 |
# Infected Departures --------------------------------------------------------- |
|
40 | 11x |
res <- update_active(rate, rand, active, status, label = "i") |
41 | 11x |
nDepartures <- res[[1]] |
42 | ! |
if(!is.null(res[[2]])) active <- dat$attr$active[res[[2]]] <- 0 |
43 |
|
|
44 | ! |
if (at == 2) dat$epi$di.flow <- c(0, nDepartures) else dat$epi$di.flow[at] <- nDepartures |
45 |
|
|
46 |
# Quarantined Departures --------------------------------------------------------- |
|
47 | 11x |
res <- update_active(rate, rand, active, status, label = "q") |
48 | 11x |
nDepartures <- res[[1]] |
49 | ! |
if(!is.null(res[[2]])) active <- dat$attr$active[res[[2]]] <- 0 |
50 |
|
|
51 | ! |
if (at == 2) dat$epi$dq.flow <- c(0, nDepartures) else dat$epi$dq.flow[at] <- nDepartures |
52 |
|
|
53 |
# Hospitalised Departures --------------------------------------------------------- |
|
54 | 11x |
res <- update_active(rate, rand, active, status, label = "h") |
55 | 11x |
nDepartures <- res[[1]] |
56 | ! |
if(!is.null(res[[2]])) active <- dat$attr$active[res[[2]]] <- 0 |
57 |
|
|
58 | ! |
if (at == 2) dat$epi$dh.flow <- c(0, nDepartures) else dat$epi$dh.flow[at] <- nDepartures |
59 |
|
|
60 |
# Recovered Departures -------------------------------------------------------- |
|
61 | 11x |
res <- update_active(rate, rand, active, status, label = "r") |
62 | 11x |
nDepartures <- res[[1]] |
63 | ! |
if(!is.null(res[[2]])) active <- dat$attr$active[res[[2]]] <- 0 |
64 |
|
|
65 | ! |
if (at == 2) dat$epi$dr.flow <- c(0, nDepartures) else dat$epi$dr.flow[at] <- nDepartures |
66 |
|
|
67 |
# return -------------------------------------------------------- |
|
68 | 11x |
return(dat) |
69 |
} |
1 |
#' @title Initial Conditions for Stochastic Individual Contact Models |
|
2 |
#' |
|
3 |
#' @description Sets the initial conditions for stochastic individual contact |
|
4 |
#' models simulated with \code{icm}. |
|
5 |
#' |
|
6 |
#' @param s.num Initial number of *S compartment individuals in |
|
7 |
#' the simulated population. An overall population of 10,000 is a good |
|
8 |
#' compromise. A set of models will still take several minutes or more |
|
9 |
#' to run, in parallel. |
|
10 |
#' @param e.num Initial number of E compartment individuals in |
|
11 |
#' the simulated population. |
|
12 |
#' @param i.num Initial number of I compartment individuals in |
|
13 |
#' the simulated population. |
|
14 |
#' @param q.num Initial number of Q compartment individuals in |
|
15 |
#' the simulated population. |
|
16 |
#' @param h.num Initial number of H compartment individuals in |
|
17 |
#' the simulated population. |
|
18 |
#' @param r.num Initial number of R compartment individuals in |
|
19 |
#' the simulated population. |
|
20 |
#' @param f.num Initial number of F compartment individuals in |
|
21 |
#' the simulated population. |
|
22 |
#' @param ... Additional initial conditions passed to model. |
|
23 |
#' |
|
24 |
#' @details |
|
25 |
#' The initial conditions for a model solved with \code{\link{icm}} should be |
|
26 |
#' input into the \code{init.icm} function. This function handles initial |
|
27 |
#' conditions for both base models and original models using new modules. For |
|
28 |
#' an overview of initial conditions for base ICM class models, consult the |
|
29 |
#' \href{http://statnet.github.io/tut/BasicICMs.html}{Basic ICMs} tutorial. |
|
30 |
#' |
|
31 |
#' @seealso Use \code{\link{param.icm}} to specify model parameters and |
|
32 |
#' \code{\link{control.icm}} to specify the control settings. Run the |
|
33 |
#' parameterized model with \code{\link{icm}}. |
|
34 |
#' |
|
35 |
#' @keywords parameterization |
|
36 |
#' |
|
37 |
#' @export |
|
38 |
#' |
|
39 | ||
40 |
init_seiqhrf <- function(s.num = 9997, e.num=0, i.num = 3, q.num=0, h.num=0, r.num = 0, f.num = 0, ...) { |
|
41 | ||
42 | 5x |
p <- list() |
43 | 5x |
passed <- names(as.list(match.call())[-1]) |
44 | 5x |
p$usr.specified <- passed |
45 |
|
|
46 | 5x |
formal.args <- formals(sys.function()) |
47 | 5x |
formal.args[["..."]] <- NULL |
48 | 5x |
for (arg in names(formal.args)) { |
49 | 35x |
if (as.logical(mget(arg) != "")) { |
50 | 35x |
p[arg] <- list(get(arg)) |
51 |
} |
|
52 |
} |
|
53 | 5x |
dot.args <- list(...) |
54 | 5x |
names.dot.args <- names(dot.args) |
55 | 5x |
if (length(dot.args) > 0) { |
56 | ! |
for (i in 1:length(dot.args)) { |
57 | ! |
p[[names.dot.args[i]]] <- dot.args[[i]] |
58 |
} |
|
59 |
} |
|
60 |
|
|
61 |
|
|
62 |
## Output |
|
63 | 5x |
class(p) <- c("init.seiqhrf", "list") |
64 | 5x |
return(p) |
65 |
} |
1 |
#' Control Settings for Stochastic Individual Contact Models |
|
2 |
#' |
|
3 |
#' Sets the controls for stochastic individual contact models simulated |
|
4 |
#' with \code{\link{simulate_seiqhrf}}. Similar to EpiModel::control.icm, but allows for |
|
5 |
#' model types with additional compartments (e.g. 'SEIQHRF'). |
|
6 |
#' |
|
7 |
#' @param type Disease type to be modeled, with the choice of \code{"SI"} for |
|
8 |
#' Susceptible-Infected diseases, \code{"SIR"} for |
|
9 |
#' Susceptible-Infected-Recovered diseases, and \code{"SIS"} for |
|
10 |
#' Susceptible-Infected-Susceptible diseases. |
|
11 |
#' @param nsteps Number of time steps to solve the model over. This must be a |
|
12 |
#' positive integer. |
|
13 |
#' @param nsims Number of simulations to run. |
|
14 |
#' @param prog.rand Method for progression from E compartment to I. If TRUE, |
|
15 |
#' random binomial draws at prog.rate, if FALSE, random draws from a |
|
16 |
#' Weibull distribution (yes, I know it should be a discrete Weibull |
|
17 |
#' distribution but it makes little difference and speed of computation |
|
18 |
#' matters), with parameters prog.dist.scale and prog.dist.shape |
|
19 |
#' @param rec.rand If \code{TRUE}, use a stochastic recovery model, with the |
|
20 |
#' number of recovered at each time step a function of random draws from |
|
21 |
#' a binomial distribution with the probability equal to \code{rec.rate}. |
|
22 |
#' If \code{FALSE}, then a deterministic rounded count of the expectation |
|
23 |
#' implied by that rate. |
|
24 |
#' @param quar.rand Method for self-isolation transition from I to Q. If TRUE, |
|
25 |
#' random binomial draws at quar.rate, if FALSE, random sample with a |
|
26 |
#' sample fraction also given by `quar.rate. |
|
27 |
#' @param hosp.rand Method for transition from I or Q to H -- that is, from |
|
28 |
#' infectious or from self-isolated to requiring hospitalisation. If |
|
29 |
#' TRUE, random binomial draws at hosp.rate, if FALSE, random sample |
|
30 |
#' with a sample fraction also given by `hosp.rate. |
|
31 |
#' @param disch.rand Method for transition from H to R -- that is, from |
|
32 |
#' requiring hospitalisation to recovered. If TRUE, random binomial |
|
33 |
#' draws at disch.rate, if FALSE, random sample with a sample fraction |
|
34 |
#' also given by disch.rate. Note that the only way out of the H |
|
35 |
#' compartment is recovery or death. |
|
36 |
#' @param fat.rand Method for case fatality transition from H to F. If TRUE, |
|
37 |
#' random binomial draws at fat.rate.base, if FALSE, random sample with |
|
38 |
#' a sample fraction also given by fat.rate.base. However, if the |
|
39 |
#' current number of patients in the H (needs hospitalisation) |
|
40 |
#' compartment is above a hospital capacity level specified by |
|
41 |
#' hosp.cap, then the fatality rate is the mean of the base fatality |
|
42 |
#' rate weighted by the hospital capacity, plus a higher rate, |
|
43 |
#' specified by fat.rate.overcap, weighted by the balance of those |
|
44 |
#' requiring hospitalisation (but who can't be accommodated). By |
|
45 |
#' setting fat.rate.overcap higher, the effect of exceeding the |
|
46 |
#' capacity of the health care system can be simulated. There is also |
|
47 |
#' a coefficient fat.tcoeff for the fatality rates that increases them |
|
48 |
#' as a linear function of the number of days the individual has been |
|
49 |
#' in the H compartment. Use of the co-efficient better approximates |
|
50 |
#' the trapezoid survival time distribution typical of ICU patients. |
|
51 |
#' @param arec.rand Method for recovery transition from E to R. If TRUE, |
|
52 |
#' random binomial draws at arec.rate, if FALSE, random draws from a |
|
53 |
#' random draws from a Weibull distribution, with parameters |
|
54 |
#' arec.dist.scale and arec.dist.shape. |
|
55 |
#' @param a.rand If \code{TRUE}, use a stochastic arrival model, with the |
|
56 |
#' number of arrivals at each time step a function of random draws from a |
|
57 |
#' binomial distribution with the probability equal to the governing arrival |
|
58 |
#' rates. If \code{FALSE}, then a deterministic rounded count of the |
|
59 |
#' expectation implied by those rates. |
|
60 |
#' @param d.rand If \code{TRUE}, use a stochastic departure model, with the number of |
|
61 |
#' departures at each time step a function of random draws from a binomial |
|
62 |
#' distribution with the probability equal to the governing departure rates. |
|
63 |
#' If \code{FALSE}, then a deterministic rounded count of the expectation |
|
64 |
#' implied by those rates. |
|
65 |
#' @param initialize.FUN Module to initialize the model at the outset, with the |
|
66 |
#' default function of \code{\link{initialize.icm}}. |
|
67 |
#' @param infection.FUN Module to simulate disease infection, with the default |
|
68 |
#' function of \code{\link{infection.icm}}. |
|
69 |
#' @param recovery.FUN Module to simulate disease recovery, with the default |
|
70 |
#' function of \code{\link{recovery.icm}}. |
|
71 |
#' @param departures.FUN Module to simulate departures or exits, with the default |
|
72 |
#' function of \code{\link{departures.icm}}. |
|
73 |
#' @param arrivals.FUN Module to simulate arrivals or entries, with the default |
|
74 |
#' function of \code{\link{arrivals.icm}}. |
|
75 |
#' @param get_prev.FUN Module to calculate disease prevalence at each time step, |
|
76 |
#' with the default function of \code{\link{get_prev.icm}}. |
|
77 |
#' @param verbose If \code{TRUE}, print model progress to the console. |
|
78 |
#' @param verbose.int Time step interval for printing progress to console, where |
|
79 |
#' 0 (the default) prints completion status of entire simulation and |
|
80 |
#' positive integer \code{x} prints progress after each \code{x} time |
|
81 |
#' steps. |
|
82 |
#' @param skip.check If \code{TRUE}, skips the default error checking for the |
|
83 |
#' structure and consistency of the parameter values, initial conditions, |
|
84 |
#' and control settings before running base epidemic models. Setting |
|
85 |
#' this to \code{FALSE} is recommended when running models with new modules |
|
86 |
#' specified. |
|
87 |
#' @param ncores Number of physical CPU cores used for parallel computation. |
|
88 |
#' @param ... Additional control settings passed to model. |
|
89 |
#' |
|
90 |
#' @details |
|
91 |
#' \code{control} sets the required control settings for any stochastic |
|
92 |
#' individual contact model solved with the \code{simulate_seiqhrf} function. Controls |
|
93 |
#' are required for both base model types and when passing original process |
|
94 |
#' modules. For an overview of control settings for base ICM class models, |
|
95 |
#' consult the \href{http://statnet.github.io/tut/BasicICMs.html}{Basic ICMs} |
|
96 |
#' tutorial. For all base models, the \code{type} argument is a necessary |
|
97 |
#' parameter and it has no default. |
|
98 |
#' |
|
99 |
#' @section New Modules: |
|
100 |
#' Base ICM models use a set of module functions that specify |
|
101 |
#' how the individual agents in the population are subjected to infection, recovery, |
|
102 |
#' demographics, and other processes. Core modules are those listed in the |
|
103 |
#' \code{.FUN} arguments. For each module, there is a default function used in |
|
104 |
#' the simulation. The default infection module, for example, is contained in |
|
105 |
#' the \code{\link{infection.FUN}} function. |
|
106 |
#' |
|
107 |
#' For original models, one may substitute replacement module functions for any of |
|
108 |
#' the default functions. New modules may be added to the workflow at each time |
|
109 |
#' step by passing a module function via the \code{...} argument. |
|
110 |
#' |
|
111 |
#' @keywords parameterization |
|
112 |
#' |
|
113 |
#' @return A list of control parameters and core functions |
|
114 |
#' @export |
|
115 |
control_seiqhrf <- function(type = "SEIQHRF", nsteps = 366, nsims = 8, |
|
116 |
prog.rand = FALSE, quar.rand = TRUE, hosp.rand = TRUE, |
|
117 |
disch.rand = TRUE, rec.rand = FALSE, arec.rand = TRUE, |
|
118 |
fat.rand = TRUE, a.rand = TRUE, d.rand = TRUE, |
|
119 |
initialize.FUN = 'initialize.FUN', infection.FUN = 'infection.FUN', |
|
120 |
recovery.FUN = 'recovery.FUN', departures.FUN = 'departures.FUN', |
|
121 |
arrivals.FUN = 'arrivals.FUN', get_prev.FUN = 'get_prev.FUN', |
|
122 |
verbose = FALSE, verbose.int = 0, skip.check = FALSE, |
|
123 |
ncores = 4, ...) { |
|
124 | ||
125 |
# Get arguments |
|
126 | 5x |
p <- list() |
127 | 5x |
passed <- names(as.list(match.call())[-1]) |
128 | 5x |
p$usr.specified <- passed |
129 |
|
|
130 | 5x |
formal.args <- formals(sys.function()) |
131 | 5x |
formal.args[["..."]] <- NULL |
132 | 5x |
for (arg in names(formal.args)) { |
133 | 110x |
if (as.logical(mget(arg) != "")) { |
134 | 110x |
p[arg] <- list(get(arg)) |
135 |
} |
|
136 |
} |
|
137 | 5x |
dot.args <- list(...) |
138 | 5x |
names.dot.args <- names(dot.args) |
139 | 5x |
if (length(dot.args) > 0) { |
140 | ! |
for (i in 1:length(dot.args)) { |
141 | ! |
p[[names.dot.args[i]]] <- dot.args[[i]] |
142 |
} |
|
143 |
} |
|
144 | ||
145 | 5x |
if ("births.FUN" %in% names(dot.args)) { |
146 | ! |
p$arrivals.FUN <- dot.args$births.FUN |
147 | ! |
p$births.FUN <- dot.args$births.FUN <- NULL |
148 | ! |
message("EpiModel 1.7.0 onward renamed the birth function births.FUN to |
149 | ! |
arrivals.FUN. See documentation for details.") |
150 |
} |
|
151 | 5x |
if ("deaths.FUN" %in% names(dot.args)) { |
152 | ! |
p$departures.FUN <- dot.args$deaths.FUN |
153 | ! |
p$deaths.FUN <- dot.args$deaths.FUN <- NULL |
154 | ! |
message("EpiModel 1.7.0 onward renamed the death function deaths.FUN to |
155 | ! |
departures.FUN. See documentation for details.") |
156 |
} |
|
157 | ||
158 | ||
159 |
## Module classification |
|
160 | 5x |
p$bi.mods <- grep(".FUN", names(formal.args), value = TRUE) |
161 |
#p$user.mods <- grep(".FUN", names(dot.args), value = TRUE) |
|
162 | ||
163 | ||
164 |
## Defaults and checks |
|
165 | 5x |
if (is.null(p$type) | !(p$type %in% c("SI", "SIS", "SIR", "SEIR", "SEIQHR", |
166 | 5x |
"SEIQHRF"))) { |
167 | ! |
stop("Specify type as \"SI\", \"SIS\", \"SIR\", \"SEIR\", \"SEIQHR\", or |
168 | ! |
\"SEIQHRF\" ", call. = FALSE) |
169 |
} |
|
170 | ! |
if (is.null(p$nsteps)) stop("Specify nsteps", call. = FALSE) |
171 | ! |
if (is.null(p$nsims)) stop("Specify nsims", call. = FALSE) |
172 | ||
173 |
## Output |
|
174 | 5x |
class(p) <- c("control.seiqhrf", "list") |
175 | 5x |
return(p) |
176 |
} |
1 |
#' Get previous function |
|
2 |
#' |
|
3 |
#' Utility function that collects prevalence and transition time data from each |
|
4 |
#' run and stores it away in the simulation result object. |
|
5 |
#' |
|
6 |
#' @param dat Merged input parameters. |
|
7 |
#' @param at Step number |
|
8 |
#' |
|
9 |
#' @return Updated dat |
|
10 |
#' @importFrom EpiModel get_prev.icm |
|
11 |
#' @export |
|
12 |
get_prev.FUN <- function(dat, at) { |
|
13 | ||
14 | ! |
if (at == 1) { |
15 | ! |
dat$epi <- list() |
16 | ! |
dat$epi$s.num <- sum(dat$attr$active == 1 & |
17 | ! |
dat$attr$status == "s" & |
18 | ! |
dat$attr$group == 1) |
19 | ! |
dat$epi$i.num <- sum(dat$attr$active == 1 & |
20 | ! |
dat$attr$status == "i" & |
21 | ! |
dat$attr$group == 1) |
22 | ! |
if (dat$control$type %in% c("SEIR", "SEIQHR", "SEIQHRF")) { |
23 | ! |
dat$epi$e.num <- sum(dat$attr$active == 1 & |
24 | ! |
dat$attr$status == "e" & |
25 | ! |
dat$attr$group == 1) |
26 |
} |
|
27 | ! |
if (dat$control$type %in% c("SIR", "SEIR", "SEIQHR", "SEIQHRF")) { |
28 | ! |
dat$epi$r.num <- sum(dat$attr$active == 1 & |
29 | ! |
dat$attr$status == "r" & |
30 | ! |
dat$attr$group == 1) |
31 |
} |
|
32 | ! |
if (dat$control$type %in% c("SEIQHR", "SEIQHRF")) { |
33 | ! |
dat$epi$q.num <- sum(dat$attr$active == 1 & |
34 | ! |
dat$attr$status == "q" & |
35 | ! |
dat$attr$group == 1) |
36 | ! |
dat$epi$h.num <- sum(dat$attr$active == 1 & |
37 | ! |
dat$attr$status == "h" & |
38 | ! |
dat$attr$group == 1) |
39 |
} |
|
40 | ! |
if (dat$control$type =="SEIQHRF") { |
41 | ! |
dat$epi$f.num <- sum(dat$attr$active == 1 & |
42 | ! |
dat$attr$status == "f" & |
43 | ! |
dat$attr$group == 1) |
44 |
} |
|
45 | ! |
if (dat$control$type == "SIR") { |
46 | ! |
dat$epi$num <- dat$epi$s.num + |
47 | ! |
dat$epi$i.num + |
48 | ! |
dat$epi$r.num |
49 | ! |
} else if (dat$control$type == "SEIR") { |
50 | ! |
dat$epi$num <- dat$epi$s.num + |
51 | ! |
dat$epi$e.num + |
52 | ! |
dat$epi$i.num + |
53 | ! |
dat$epi$r.num |
54 | ! |
} else if (dat$control$type == "SEIQHR") { |
55 | ! |
dat$epi$num <- dat$epi$s.num + |
56 | ! |
dat$epi$e.num + |
57 | ! |
dat$epi$i.num + |
58 | ! |
dat$epi$q.num + |
59 | ! |
dat$epi$h.num + |
60 | ! |
dat$epi$r.num |
61 | ! |
} else if (dat$control$type == "SEIQHRF") { |
62 | ! |
dat$epi$num <- dat$epi$s.num + |
63 | ! |
dat$epi$e.num + |
64 | ! |
dat$epi$i.num + |
65 | ! |
dat$epi$q.num + |
66 | ! |
dat$epi$h.num + |
67 | ! |
dat$epi$r.num + |
68 | ! |
dat$epi$f.num |
69 |
} else { |
|
70 | ! |
dat$epi$num <- dat$epi$s.num + dat$epi$i.num |
71 |
} |
|
72 | ! |
if (dat$param$groups == 2) { |
73 | ! |
dat$epi$s.num.g2 <- sum(dat$attr$active == 1 & |
74 | ! |
dat$attr$status == "s" & |
75 | ! |
dat$attr$group == 2) |
76 | ! |
if (dat$control$type %in% c("SEIR", "SEIQHR", "SEIQHRF")) { |
77 | ! |
dat$epi$e.num.g2 <- sum(dat$attr$active == 1 & |
78 | ! |
dat$attr$status == "e" & |
79 | ! |
dat$attr$group == 2) |
80 |
} |
|
81 | ! |
if (dat$control$type %in% c("SEIQHR", "SEIQHRF")) { |
82 | ! |
dat$epi$q.num.g2 <- sum(dat$attr$active == 1 & |
83 | ! |
dat$attr$status == "q" & |
84 | ! |
dat$attr$group == 2) |
85 | ! |
dat$epi$h.num.g2 <- sum(dat$attr$active == 1 & |
86 | ! |
dat$attr$status == "h" & |
87 | ! |
dat$attr$group == 2) |
88 |
} |
|
89 | ! |
if (dat$control$type %in% c("SEIQHRF")) { |
90 | ! |
dat$epi$f.num.g2 <- sum(dat$attr$active == 1 & |
91 | ! |
dat$attr$status == "f" & |
92 | ! |
dat$attr$group == 2) |
93 |
} |
|
94 | ! |
dat$epi$i.num.g2 <- sum(dat$attr$active == 1 & |
95 | ! |
dat$attr$status == "i" & |
96 | ! |
dat$attr$group == 2) |
97 | ! |
dat$epi$num.g2 <- dat$epi$s.num.g2 + dat$epi$i.num.g2 |
98 | ! |
if (dat$control$type %in% c("SIR", "SEIR", "SEIQHR", "SEIQHRF")) { |
99 | ! |
dat$epi$r.num.g2 <- sum(dat$attr$active == 1 & |
100 | ! |
dat$attr$status == "r" & |
101 | ! |
dat$attr$group == 2) |
102 |
} |
|
103 | ! |
if (dat$control$type == "SIR") { |
104 | ! |
dat$epi$num.g2 <- dat$epi$s.num.g2 + |
105 | ! |
dat$epi$i.num.g2 + |
106 | ! |
dat$epi$r.num.g2 |
107 | ! |
} else if (dat$control$type == "SEIR") { |
108 | ! |
dat$epi$num.g2 <- dat$epi$s.num.g2 + |
109 | ! |
dat$epi$e.num.g2 + |
110 | ! |
dat$epi$i.num.g2 + |
111 | ! |
dat$epi$r.num.g2 |
112 | ! |
} else if (dat$control$type == "SEIQHR") { |
113 | ! |
dat$epi$num.g2 <- dat$epi$s.num.g2 + |
114 | ! |
dat$epi$e.num.g2 + |
115 | ! |
dat$epi$i.num.g2 + |
116 | ! |
dat$epi$q.num.g2 + |
117 | ! |
dat$epi$h.num.g2 + |
118 | ! |
dat$epi$r.num.g2 |
119 | ! |
} else if (dat$control$type == "SEIQHRF") { |
120 | ! |
dat$epi$num.g2 <- dat$epi$s.num.g2 + |
121 | ! |
dat$epi$e.num.g2 + |
122 | ! |
dat$epi$i.num.g2 + |
123 | ! |
dat$epi$q.num.g2 + |
124 | ! |
dat$epi$h.num.g2 + |
125 | ! |
dat$epi$r.num.g2 + |
126 | ! |
dat$epi$f.num.g2 |
127 |
} else { |
|
128 | ! |
dat$epi$num.g2 <- dat$epi$s.num.g2 + dat$epi$i.num.g2 |
129 |
} |
|
130 |
} |
|
131 |
} else { |
|
132 | ! |
dat$epi$s.num[at] <- sum(dat$attr$active == 1 & |
133 | ! |
dat$attr$status == "s" & |
134 | ! |
dat$attr$group == 1) |
135 | ! |
if (dat$control$type %in% c("SEIR", "SEIQHR", "SEIQHRF")) { |
136 | ! |
dat$epi$e.num[at] <- sum(dat$attr$active == 1 & |
137 | ! |
dat$attr$status == "e" & |
138 | ! |
dat$attr$group == 1) |
139 |
} |
|
140 | ! |
dat$epi$i.num[at] <- sum(dat$attr$active == 1 & |
141 | ! |
dat$attr$status == "i" & |
142 | ! |
dat$attr$group == 1) |
143 | ! |
if (dat$control$type %in% c("SIR", "SEIR", "SEIQHR", "SEIQHRF")) { |
144 | ! |
dat$epi$r.num[at] <- sum(dat$attr$active == 1 & |
145 | ! |
dat$attr$status == "r" & |
146 | ! |
dat$attr$group == 1) |
147 |
} |
|
148 | ! |
if (dat$control$type %in% c("SEIQHR", "SEIQHRF")) { |
149 | ! |
dat$epi$q.num[at] <- sum(dat$attr$active == 1 & |
150 | ! |
dat$attr$status == "q" & |
151 | ! |
dat$attr$group == 1) |
152 | ! |
dat$epi$h.num[at] <- sum(dat$attr$active == 1 & |
153 | ! |
dat$attr$status == "h" & |
154 | ! |
dat$attr$group == 1) |
155 |
} |
|
156 | ! |
if (dat$control$type %in% c("SEIQHRF")) { |
157 | ! |
dat$epi$f.num[at] <- sum(dat$attr$active == 1 & |
158 | ! |
dat$attr$status == "f" & |
159 | ! |
dat$attr$group == 1) |
160 |
} |
|
161 | ! |
if (dat$control$type == "SIR") { |
162 | ! |
dat$epi$num[at] <- dat$epi$s.num[at] + |
163 | ! |
dat$epi$i.num[at] + |
164 | ! |
dat$epi$r.num[at] |
165 | ! |
} else if (dat$control$type == "SEIR") { |
166 | ! |
dat$epi$num[at] <- dat$epi$s.num[at] + |
167 | ! |
dat$epi$e.num[at] + |
168 | ! |
dat$epi$i.num[at] + |
169 | ! |
dat$epi$r.num[at] |
170 | ! |
} else if (dat$control$type == "SEIQHR") { |
171 | ! |
dat$epi$num[at] <- dat$epi$s.num[at] + |
172 | ! |
dat$epi$e.num[at] + |
173 | ! |
dat$epi$i.num[at] + |
174 | ! |
dat$epi$q.num[at] + |
175 | ! |
dat$epi$h.num[at] + |
176 | ! |
dat$epi$r.num[at] |
177 | ! |
} else if (dat$control$type == "SEIQHRF") { |
178 | ! |
dat$epi$num[at] <- dat$epi$s.num[at] + |
179 | ! |
dat$epi$e.num[at] + |
180 | ! |
dat$epi$i.num[at] + |
181 | ! |
dat$epi$q.num[at] + |
182 | ! |
dat$epi$h.num[at] + |
183 | ! |
dat$epi$r.num[at] + |
184 | ! |
dat$epi$f.num[at] |
185 |
} else { |
|
186 | ! |
dat$epi$num[at] <- dat$epi$s.num[at] + dat$epi$i.num[at] |
187 |
} |
|
188 | ||
189 | ! |
if (dat$param$groups == 2) { |
190 | ! |
dat$epi$s.num.g2[at] <- sum(dat$attr$active == 1 & |
191 | ! |
dat$attr$status == "s" & |
192 | ! |
dat$attr$group == 2) |
193 | ! |
if (dat$control$type %in% c("SEIR", "SEIQHR", "SEIQHRF")) { |
194 | ! |
dat$epi$i.num.g2[at] <- sum(dat$attr$active == 1 & |
195 | ! |
dat$attr$status == "e" & |
196 | ! |
dat$attr$group == 2) |
197 |
} |
|
198 | ! |
dat$epi$i.num.g2[at] <- sum(dat$attr$active == 1 & |
199 | ! |
dat$attr$status == "i" & |
200 | ! |
dat$attr$group == 2) |
201 | ! |
if (dat$control$type %in% c("SIR", "SEIR", "SEIQHR", "SEIQHRF")) { |
202 | ! |
dat$epi$r.num.g2[at] <- sum(dat$attr$active == 1 & |
203 | ! |
dat$attr$status == "r" & |
204 | ! |
dat$attr$group == 2) |
205 |
} |
|
206 | ! |
if (dat$control$type %in% c("SEIQHR", "SEIQHRF")) { |
207 | ! |
dat$epi$q.num.g2[at] <- sum(dat$attr$active == 1 & |
208 | ! |
dat$attr$status == "q" & |
209 | ! |
dat$attr$group == 2) |
210 | ! |
dat$epi$h.num.g2[at] <- sum(dat$attr$active == 1 & |
211 | ! |
dat$attr$status == "h" & |
212 | ! |
dat$attr$group == 2) |
213 |
} |
|
214 | ! |
if (dat$control$type %in% c("SEIQHRF")) { |
215 | ! |
dat$epi$f.num.g2[at] <- sum(dat$attr$active == 1 & |
216 | ! |
dat$attr$status == "f" & |
217 | ! |
dat$attr$group == 2) |
218 |
} |
|
219 | ! |
if (dat$control$type == "SIR") { |
220 | ! |
dat$epi$num.g2[at] <- dat$epi$s.num.g2[at] + |
221 | ! |
dat$epi$i.num.g2[at] + |
222 | ! |
dat$epi$r.num.g2[at] |
223 | ! |
} else if (dat$control$type == "SEIR") { |
224 | ! |
dat$epi$num.g2[at] <- dat$epi$s.num.g2[at] + |
225 | ! |
dat$epi$e.num.g2[at] + |
226 | ! |
dat$epi$i.num.g2[at] + |
227 | ! |
dat$epi$r.num.g2[at] |
228 | ! |
} else if (dat$control$type == "SEIQHR") { |
229 | ! |
dat$epi$num.g2[at] <- dat$epi$s.num.g2[at] + |
230 | ! |
dat$epi$e.num.g2[at] + |
231 | ! |
dat$epi$i.num.g2[at] + |
232 | ! |
dat$epi$q.num.g2[at] + |
233 | ! |
dat$epi$h.num.g2[at] + |
234 | ! |
dat$epi$r.num.g2[at] |
235 | ! |
} else if (dat$control$type == "SEIQHRF") { |
236 | ! |
dat$epi$num.g2[at] <- dat$epi$s.num.g2[at] + |
237 | ! |
dat$epi$e.num.g2[at] + |
238 | ! |
dat$epi$i.num.g2[at] + |
239 | ! |
dat$epi$q.num.g2[at] + |
240 | ! |
dat$epi$h.num.g2[at] + |
241 | ! |
dat$epi$r.num.g2[at] + |
242 | ! |
dat$epi$f.num.g2[at] |
243 |
} else { |
|
244 | ! |
dat$epi$num.g2[at] <- dat$epi$s.num.g2[at] + |
245 | ! |
dat$epi$i.num.g2[at] |
246 |
} |
|
247 |
} |
|
248 |
} |
|
249 | ||
250 | ! |
return(dat) |
251 |
} |
1 |
#' @title Epidemic Parameters for Stochastic Individual Contact Models |
|
2 |
#' |
|
3 |
#' @description Sets the epidemic parameters for stochastic individual contact |
|
4 |
#' models simulated with \code{seiqhrf}. |
|
5 |
#' @param inf.prob.e Probability of passing on infection at each |
|
6 |
#' exposure event for interactions between infectious people in the E |
|
7 |
#' compartment and susceptible in S. Note the default is lower than for |
|
8 |
#' inf.prob.i reflecting the reduced infectivity of infected but |
|
9 |
#' asymptomatic people (the E compartment). Otherwise as for inf.exp.i. |
|
10 |
#' @param act.rate.e The number of exposure events (acts) between |
|
11 |
#' infectious individuals in the E compartment and susceptible individuals |
|
12 |
#' in the S compartment, per day. Otherwise as for act.rate.i. |
|
13 |
#' @param inf.prob.i Probability of passing on infection at each |
|
14 |
#' exposure event for interactions between infectious people in the I |
|
15 |
#' compartment and susceptible in S. Reducing inf.prob.i is equivalent to |
|
16 |
#' increasing hygiene measures, such as not putting hands in eyes, nose or |
|
17 |
#' moth, use of hand sanitisers, wearing masks by the infected, and so on. |
|
18 |
#' @param act.rate.i The number of exposure events (acts) between |
|
19 |
#' infectious individuals in the I compartment and susceptible individuals |
|
20 |
#' in the S compartment, per day. It's stochastic, so the rate is an |
|
21 |
#' average, some individuals may have more or less. Note that not every |
|
22 |
#' exposure event results in infection - that is governed by the inf.prob.i |
|
23 |
#' parameters (see below). Reducing act.rate.i is equivalent to increasing |
|
24 |
#' social distancing by people in the I compartment. |
|
25 |
#' @param inf.prob.q Probability of passing on infection at each |
|
26 |
#' exposure event for interactions between infectious people in the Q |
|
27 |
#' compartment and susceptible in S. Note the default is lower than for |
|
28 |
#' inf.prob.i reflecting the greater care that self-isolated individuals |
|
29 |
#' will, on average, take regarding hygiene measures, such as wearing masks, |
|
30 |
#' to limit spread to others. Otherwise as for inf.exp.i. |
|
31 |
#' @param act.rate.q The number of exposure events (acts) between |
|
32 |
#' infectious individuals in the Q compartment (isolated, self or otherwise) |
|
33 |
#' and susceptible individuals in the S compartment, per day. Note the much |
|
34 |
#' lower rate than for the I and E compartments, reflecting the much |
|
35 |
#' greater degree of social isolation for someone in (self-)isolation. The |
|
36 |
#' exposure event rate is not zero for this group, just much less. |
|
37 |
#' Otherwise as for act.rate.i. |
|
38 |
#' @param quar.rate Rate per day at which symptomatic (or tested |
|
39 |
#' positive), infected I compartment people enter self-isolation (Q |
|
40 |
#' compartment). Asymptomatic E compartment people can't enter |
|
41 |
#' self-isolation because they don't yet know they are infected. Default is |
|
42 |
#' a low rate reflecting low community awareness or compliance with |
|
43 |
#' self-isolation requirements or practices, but this can be tweaked when |
|
44 |
#' exploring scenarios. |
|
45 |
#' @param quar.dist.scale Scale parameter for Weibull distribution for |
|
46 |
#' recovery, see quar.rand for details. |
|
47 |
#' @param quar.dist.shape Shape parameter for Weibull distribution for |
|
48 |
#' recovery, see quar.rand for details. Read up on the Weibull distribution |
|
49 |
#' before changing the default. |
|
50 |
#' @param hosp.rate Rate per day at which symptomatic (or tested |
|
51 |
#' positive), infected I compartment people or self-isolated Q compartment |
|
52 |
#' people enter the state of requiring hospital care -- that is, become |
|
53 |
#' serious cases. A default rate of 1% per day with an average illness |
|
54 |
#' duration of about 10 days means a bit less than 10% of cases will |
|
55 |
#' require hospitalisation, which seems about right (but can be tweaked, |
|
56 |
#' of course). |
|
57 |
#' @param hosp.dist.scale Scale parameter for Weibull distribution for |
|
58 |
#' recovery, see hosp.rand for details. |
|
59 |
#' @param hosp.dist.shape Shape parameter for Weibull distribution for |
|
60 |
#' recovery, see hosp.rand for details. Read up on the Weibull distribution |
|
61 |
#' before changing the default. |
|
62 |
#' @param disch.rate Rate per day at which people needing |
|
63 |
#' hospitalisation recover. |
|
64 |
#' @param disch.dist.scale Scale parameter for Weibull distribution for |
|
65 |
#' recovery, see disch.rand for details. |
|
66 |
#' @param disch.dist.shape Shape parameter for Weibull distribution for |
|
67 |
#' recovery, see disch.rand for details. Read up on the Weibull distribution |
|
68 |
#' before changing the default. |
|
69 |
#' @param prog.rate Rate per day at which people who are infected |
|
70 |
#' but asymptomatic (E compartment) progress to becoming symptomatic (or |
|
71 |
#' test-positive), the I compartment. See prog.rand above for more details. |
|
72 |
#' @param prog.dist.scale Scale parameter for Weibull distribution |
|
73 |
#' for progression, see prog.rand for details. |
|
74 |
#' @param prog.dist.shape Shape parameter for Weibull distribution |
|
75 |
#' for progression, see prog.rand for details. Read up on the Weibull |
|
76 |
#' distribution before changing the default. |
|
77 |
#' @param rec.rate Rate per day at which people who are infected and |
|
78 |
#' symptomatic (I compartment) recover, thus entering the R compartment. |
|
79 |
#' See rec.rand above for more details. |
|
80 |
#' @param rec.dist.scale Scale parameter for Weibull distribution for |
|
81 |
#' recovery, see rec.rand for details. |
|
82 |
#' @param rec.dist.shape Shape parameter for Weibull distribution for |
|
83 |
#' recovery, see rec.rand for details. Read up on the Weibull distribution |
|
84 |
#' before changing the default. |
|
85 |
#' @param arec.rate Rate per day at which people who are exposed but asymptotic |
|
86 |
#' (E compartment) recover, thus entering the R compartment. |
|
87 |
#' See arec.rand above for more details. |
|
88 |
#' @param arec.dist.scale Scale parameter for Weibull distribution for |
|
89 |
#' recovery, see arec.rand for details. |
|
90 |
#' @param arec.dist.shape Shape parameter for Weibull distribution for |
|
91 |
#' recovery, see arec.rand for details. |
|
92 |
#' @param fat.rate.base Baseline mortality rate per day for people |
|
93 |
#' needing hospitalisation (deaths due to the virus). See fat.rand for more |
|
94 |
#' details. |
|
95 |
#' @param hosp.cap Number of available hospital beds for the modelled |
|
96 |
#' population. See fat.rand for more details. |
|
97 |
#' @param fat.rate.overcap Mortality rate per day for people needing |
|
98 |
#' hospitalisation but who can't get into hospital due to the hospitals |
|
99 |
#' being full (see hosp.cap and fat.rand). The default rate is twice that |
|
100 |
#' for those who do get into hospital. |
|
101 |
#' @param fat.tcoeff Time co-efficient for increasing mortality rate |
|
102 |
#' as time in the H compartment increases for each individual in it. See |
|
103 |
#' fat.rand for details. |
|
104 |
#' @param vital Enables demographics, that is, arrivals and |
|
105 |
#' departures, to and from the simulated population. |
|
106 |
#' @param flare.inf.point (Either a vector or a scalar) Time points where there is a sudden arrival of I's. |
|
107 |
#' @param flare.inf.num (same dimension as flare.inf.point) The number of I's arriving at the specified time points in flare.inf.point. |
|
108 |
#' @param a.rate Background demographic arrival rate. Currently all |
|
109 |
#' arrivals go into the S compartment, the default is approximately the |
|
110 |
#' daily birth rate for Australia. Will be extended to cover immigration in |
|
111 |
#' future versions. |
|
112 |
#' @param a.prop.e Arrivals proportion that goes to E (immigration). |
|
113 |
#' @param a.prop.i Arrivals proportion that goes to I (immigration). |
|
114 |
#' @param a.prop.q Arrivals proportion that goes to Q (immigration). |
|
115 |
#' @param ds.rate Background demographic departure (death not due to |
|
116 |
#' virus) rates. Defaults based on Australian crude death rates. Can be |
|
117 |
#' used to model emigration as well as deaths. |
|
118 |
#' @param de.rate Background demographic departure (death not due to |
|
119 |
#' virus) rates. Defaults based on Australian crude death rates. Can be |
|
120 |
#' used to model emigration as well as deaths. |
|
121 |
#' @param di.rate Background demographic departure (death not due to |
|
122 |
#' virus) rates. Defaults based on Australian crude death rates. Can be used |
|
123 |
#' to model emigration as well as deaths. |
|
124 |
#' @param dq.rate Background demographic departure (death not due to |
|
125 |
#' virus) rates. Defaults based on Australian crude death rates. Can be used |
|
126 |
#' to model emigration as well as deaths. |
|
127 |
#' @param dh.rate Background demographic departure (death not due to |
|
128 |
#' virus) rates. Defaults based on Australian crude death rates. Can be used |
|
129 |
#' to model emigration as well as deaths. |
|
130 |
#' @param dr.rate Background demographic departure (death not due to |
|
131 |
#' virus) rates. Defaults based on Australian crude death rates. Can be used |
|
132 |
#' to model emigration as well as deaths. |
|
133 |
#' @param ... Other parameters. |
|
134 |
#' |
|
135 |
#' |
|
136 |
#' @details |
|
137 |
#' \code{param_seiqhrf} sets the epidemic parameters for the stochastic individual |
|
138 |
#' contact models simulated with the \code{\link{seiqhrf}} function. Models |
|
139 |
#' may use the base types, for which these parameters are used, or new process |
|
140 |
#' modules which may use these parameters (but not necessarily). A detailed |
|
141 |
#' description of ICM parameterization for base models is found in the |
|
142 |
#' \href{http://statnet.github.io/tut/BasicICMs.html}{Basic ICMs} tutorial. |
|
143 |
#' |
|
144 |
#' For base models, the model specification will be chosen as a result of |
|
145 |
#' the model parameters entered here and the control settings in |
|
146 |
#' \code{\link{control_seiqhrf}}. |
|
147 |
#' |
|
148 |
#' @section New Modules: |
|
149 |
#' To build original models outside of the base models, new process modules |
|
150 |
#' may be constructed to replace the existing modules or to supplement the existing |
|
151 |
#' set. These are passed into the control settings in \code{\link{control_seiqhrf}}. |
|
152 |
#' New modules may use either the existing model parameters named here, an |
|
153 |
#' original set of parameters, or a combination of both. The \code{...} allows |
|
154 |
#' the user to pass an arbitrary set of original model parameters into |
|
155 |
#' \code{param_seiqhrf}. Whereas there are strict checks with default modules for |
|
156 |
#' parameter validity, these checks are the user's responsibility with new modules. |
|
157 |
#' |
|
158 |
#' @seealso Use \code{\link{init.icm}} to specify the initial conditions and |
|
159 |
#' \code{\link{control_seiqhrf}} to specify the control settings. Run the |
|
160 |
#' parameterized model with \code{\link{seiqhrf}}. |
|
161 |
#' |
|
162 |
#' @keywords parameterization |
|
163 |
#' |
|
164 |
#' @export |
|
165 |
#' |
|
166 |
param_seiqhrf <- function(inf.prob.e = 0.02, |
|
167 |
act.rate.e = 10, |
|
168 |
inf.prob.i = 0.05, |
|
169 |
act.rate.i = 10, |
|
170 |
inf.prob.q = 0.02, |
|
171 |
act.rate.q = 2.5, |
|
172 |
prog.rate = 1/10, |
|
173 |
quar.rate = 1/30, |
|
174 |
hosp.rate = 1/100, |
|
175 |
disch.rate = 1/15, |
|
176 |
rec.rate = 0.071, |
|
177 |
arec.rate = 0.05, |
|
178 |
prog.dist.scale = 5, |
|
179 |
prog.dist.shape = 1.5, |
|
180 |
quar.dist.scale = 1, |
|
181 |
quar.dist.shape = 1, |
|
182 |
hosp.dist.scale = 1, |
|
183 |
hosp.dist.shape = 1, |
|
184 |
disch.dist.scale = 1, |
|
185 |
disch.dist.shape = 1, |
|
186 |
rec.dist.scale = 35, |
|
187 |
rec.dist.shape = 1.5, |
|
188 |
arec.dist.scale = 35, |
|
189 |
arec.dist.shape = 1.5, |
|
190 |
fat.rate.base = 1/50, |
|
191 |
hosp.cap = 40, |
|
192 |
fat.rate.overcap = 1/25, |
|
193 |
fat.tcoeff = 0.5, |
|
194 |
vital = TRUE, |
|
195 |
flare.inf.point = NULL, |
|
196 |
flare.inf.num = NULL, |
|
197 |
a.rate = (10.5/365)/1000, |
|
198 |
a.prop.e = 0.01, |
|
199 |
a.prop.i = 0.001, |
|
200 |
a.prop.q = 0.01, |
|
201 |
ds.rate = (7/365)/1000, |
|
202 |
de.rate = (7/365)/1000, |
|
203 |
di.rate = (7/365)/1000, |
|
204 |
dq.rate = (7/365)/1000, |
|
205 |
dh.rate = (20/365)/1000, |
|
206 |
dr.rate = (7/365)/1000, ...) { |
|
207 |
|
|
208 |
# Get arguments |
|
209 | 6x |
p <- list() |
210 | 6x |
passed <- names(as.list(match.call())[-1]) |
211 | 6x |
p$usr.specified <- passed |
212 |
|
|
213 | 6x |
formal.args <- formals(sys.function()) |
214 | 6x |
formal.args[["..."]] <- NULL |
215 | 6x |
for (arg in names(formal.args)) { |
216 | 246x |
if (as.logical(mget(arg) != "")) { |
217 | 246x |
p[arg] <- list(get(arg)) |
218 |
} |
|
219 |
} |
|
220 | 6x |
dot.args <- list(...) |
221 | 6x |
names.dot.args <- names(dot.args) |
222 | 6x |
if (length(dot.args) > 0) { |
223 | ! |
for (i in 1:length(dot.args)) { |
224 | ! |
p[[names.dot.args[i]]] <- dot.args[[i]] |
225 |
} |
|
226 |
} |
|
227 |
|
|
228 | 6x |
if ("b.rate" %in% names.dot.args) { |
229 | ! |
p$a.rate <- dot.args$b.rate |
230 | ! |
message("EpiModel 1.7.0 onward renamed the birth rate parameter b.rate to a.rate. ", |
231 | ! |
"See documentation for details.") |
232 |
} |
|
233 |
|
|
234 | 6x |
if ("b.rand" %in% names.dot.args) { |
235 | ! |
p$a.rand <- dot.args$b.rand |
236 | ! |
message("EpiModel 1.7.0 onward renamed the stochastic birth flag b.rand to a.rand. ", |
237 | ! |
"See documentation for details.") |
238 |
} |
|
239 |
|
|
240 |
## Defaults and checks |
|
241 | 6x |
if (is.null(p$act.rate)) { |
242 | 6x |
p$act.rate <- 1 |
243 |
} |
|
244 | 6x |
p$vital <- ifelse(!is.null(p$a.rate) | !is.null(p$ds.rate) | |
245 | 6x |
!is.null(p$di.rate) | !is.null(p$dr.rate), TRUE, FALSE) |
246 |
|
|
247 | 6x |
p$groups <- ifelse(any(grepl(".g2", names(p))) == TRUE, 2, 1) |
248 |
|
|
249 | 6x |
if (!is.null(p$inter.eff) && is.null(p$inter.start)) { |
250 | ! |
p$inter.start <- 1 |
251 |
} |
|
252 |
|
|
253 |
## Output |
|
254 | 6x |
class(p) <- c("param.seiqhrf", "list") |
255 | 6x |
return(p) |
256 |
} |
|
257 | ||
258 |
1 |
#' @title Extract Model Data for Stochastic Models |
|
2 |
#' |
|
3 |
#' @description This function extracts model simulations for objects of classes |
|
4 |
#' \code{seiqhrf} into a data frame using |
|
5 |
#' the generic \code{as.data.frame} function. |
|
6 |
#' |
|
7 |
#' @param x An \code{EpiModel} object of class \code{icm} or \code{netsim}. |
|
8 |
#' @param out Data output to data frame: \code{"mean"} for row means across |
|
9 |
#' simulations, \code{"sd"} for row standard deviations across simulations, |
|
10 |
#' \code{"qnt"} for row quantiles at the level specified in \code{qval}, |
|
11 |
#' or \code{"vals"} for values from individual simulations. |
|
12 |
#' @param sim If \code{out="vals"}, the simulation number to output. If not |
|
13 |
#' specified, then data from all simulations will be output. |
|
14 |
#' @param qval Quantile value required when \code{out="qnt"}. |
|
15 |
#' @param row.names See \code{\link{as.data.frame.default}}. |
|
16 |
#' @param optional See \code{\link{as.data.frame.default}}. |
|
17 |
#' @param ... See \code{\link{as.data.frame.default}}. |
|
18 |
#' |
|
19 |
#' @details |
|
20 |
#' These methods work for both \code{seiqhrf} class models. The |
|
21 |
#' available output includes time-specific means, standard deviations, quantiles, |
|
22 |
#' and simulation values (compartment and flow sizes) from these stochastic model |
|
23 |
#' classes. Means, standard deviations, and quantiles are calculated by taking the |
|
24 |
#' row summary (i.e., each row of data is corresponds to a time step) across all |
|
25 |
#' simulations in the model output. |
|
26 |
#' |
|
27 |
#' @method as.data.frame seiqhrf |
|
28 |
#' @keywords extract |
|
29 |
#' @importFrom stats quantile |
|
30 |
#' @export |
|
31 |
#' |
|
32 |
#' @examples |
|
33 |
#' ## Stochastic SEIQHRF model |
|
34 |
#' param <- param_seiqhrf() |
|
35 |
#' init <- init_seiqhrf(s.num = 500, i.num = 1) |
|
36 |
#' control <- control_seiqhrf(nsteps = 10, nsims = 3) |
|
37 |
#' mod <- seiqhrf(init, control, param) |
|
38 |
#' |
|
39 |
#' # Default output all simulation runs, default to all in stacked data.frame |
|
40 |
#' as.data.frame(mod) |
|
41 |
#' as.data.frame(mod, sim = 2) |
|
42 |
#' |
|
43 |
#' # Time-specific means across simulations |
|
44 |
#' as.data.frame(mod, out = "mean") |
|
45 |
#' |
|
46 |
#' # Time-specific standard deviations across simulations |
|
47 |
#' as.data.frame(mod, out = "sd") |
|
48 |
#' |
|
49 |
#' # Time-specific quantile values across simulations |
|
50 |
#' as.data.frame(mod, out = "qnt", qval = 0.25) |
|
51 |
#' as.data.frame(mod, out = "qnt", qval = 0.75) |
|
52 |
#' |
|
53 | ||
54 |
as.data.frame.seiqhrf <- function(x, row.names = NULL, optional = FALSE, |
|
55 |
out = "vals", sim, qval, ...) { |
|
56 |
|
|
57 | 10x |
df <- data.frame(time = 1:x$control$nsteps) |
58 | 10x |
nsims <- x$control$nsims |
59 |
|
|
60 | 10x |
if (out == "vals") { |
61 |
|
|
62 | ! |
if (missing(sim)) { |
63 | ! |
sim <- 1:nsims |
64 |
} |
|
65 | ! |
if (max(sim) > nsims) { |
66 | ! |
stop("Maximum sim is ", nsims, call. = FALSE) |
67 |
} |
|
68 |
|
|
69 | ! |
for (j in sim) { |
70 | ! |
if (j == min(sim)) { |
71 | ! |
for (i in seq_along(x$epi)) { |
72 | ! |
if (nsims == 1) { |
73 | ! |
df[, i + 1] <- x$epi[[i]] |
74 |
} else { |
|
75 | ! |
df[, i + 1] <- x$epi[[i]][, j] |
76 |
} |
|
77 |
} |
|
78 | ! |
df$sim <- j |
79 |
} else { |
|
80 | ! |
tdf <- data.frame(time = 1:x$control$nsteps) |
81 | ! |
for (i in seq_along(x$epi)) { |
82 | ! |
if (nsims == 1) { |
83 | ! |
tdf[, i + 1] <- x$epi[[i]] |
84 |
} else { |
|
85 | ! |
tdf[, i + 1] <- x$epi[[i]][, j] |
86 |
} |
|
87 |
} |
|
88 | ! |
tdf$sim <- j |
89 | ! |
df <- rbind(df, tdf) |
90 |
} |
|
91 |
} |
|
92 | ! |
df <- df[, c(ncol(df), 1:(ncol(df) - 1))] |
93 |
|
|
94 |
} |
|
95 |
|
|
96 |
## Output means |
|
97 | 10x |
if (out == "mean") { |
98 | 4x |
if (nsims == 1) { |
99 | ! |
for (i in seq_along(x$epi)) { |
100 | ! |
df[, i + 1] <- x$epi[[i]] |
101 |
} |
|
102 |
} |
|
103 | 4x |
if (nsims > 1) { |
104 | 4x |
for (i in seq_along(x$epi)) { |
105 | 92x |
df[, i + 1] <- rowMeans(x$epi[[i]], na.rm = TRUE) |
106 |
} |
|
107 |
} |
|
108 |
} |
|
109 |
|
|
110 |
## Output standard deviations |
|
111 | 10x |
if (out == "sd") { |
112 | 2x |
if (nsims == 1) { |
113 | ! |
for (i in seq_along(x$epi)) { |
114 | ! |
df[, i + 1] <- 0 |
115 |
} |
|
116 |
} |
|
117 | 2x |
if (nsims > 1) { |
118 | 2x |
for (i in seq_along(x$epi)) { |
119 | 46x |
df[, i + 1] <- apply(x$epi[[i]], 1, sd, na.rm = TRUE) |
120 |
} |
|
121 |
} |
|
122 |
} |
|
123 |
|
|
124 | 10x |
if (out == "qnt") { |
125 | 4x |
if (missing(qval) || length(qval) > 1 || (qval > 1 | qval < 0)) { |
126 | ! |
stop("Must specify qval as single value between 0 and 1", call. = FALSE) |
127 |
} |
|
128 | 4x |
for (i in seq_along(x$epi)) { |
129 | 92x |
df[, i + 1] <- apply(x$epi[[i]], 1, stats::quantile, probs = qval, |
130 | 92x |
na.rm = TRUE, names = FALSE) |
131 |
} |
|
132 |
} |
|
133 |
|
|
134 | 10x |
if (out == "vals") { |
135 | ! |
names(df)[3:ncol(df)] <- names(x$epi) |
136 |
} else { |
|
137 | 10x |
names(df)[2:ncol(df)] <- names(x$epi) |
138 |
} |
|
139 |
|
|
140 | 10x |
return(df) |
141 |
} |
1 | ||
2 |
#' Save icm |
|
3 |
#' |
|
4 |
#' Function to organize simulation outputs |
|
5 |
#' |
|
6 |
#' @param dat Merged input parameters. |
|
7 |
#' @param s Step number or time point. |
|
8 |
#' @param out initial out. Default = NULL |
|
9 |
#' |
|
10 |
#' @return list of output in the same format as \code{icm} object in \code{EpiModel}. |
|
11 | ||
12 |
saveout_seiqhrf <- function(dat, s, out = NULL) { |
|
13 | ||
14 | ! |
alist2df <- function(dat,s) { |
15 | ! |
alist <- list() |
16 | ! |
alist$expTime <- dat$attr$expTime |
17 | ! |
alist$infTime <- dat$attr$infTime |
18 | ! |
alist$quarTime <- dat$attr$quarTime |
19 | ! |
alist$recovTime <- dat$attr$recovTime |
20 | ! |
alist$hospTime <- dat$attr$hospTime |
21 | ! |
alist$dischTime <- dat$attr$dischTime |
22 | ! |
alist$fatTime <- dat$attr$fatTime |
23 | ! |
alist <- lapply(alist, `length<-`, max(lengths(alist))) |
24 | ! |
return(data.frame(alist)) |
25 |
} |
|
26 | ||
27 | ! |
if (s == 1) { |
28 | ! |
out <- list() |
29 | ! |
out$param <- dat$param |
30 | ! |
out$control <- dat$control |
31 | ! |
out$epi <- list() |
32 | ! |
for (j in 1:length(dat$epi)) { |
33 | ! |
out$epi[[names(dat$epi)[j]]] <- data.frame(dat$epi[j]) |
34 |
} |
|
35 |
} else { |
|
36 | ! |
for (j in 1:length(dat$epi)) { |
37 | ! |
out$epi[[names(dat$epi)[j]]][, s] <- data.frame(dat$epi[j]) |
38 |
} |
|
39 |
} |
|
40 | ||
41 |
# capture transition times from attribs |
|
42 | ! |
if (dat$control$type %in% c("SEIQHR", "SEIQHRF")) { |
43 | ! |
out$times[[paste("sim",s,sep="")]] <- alist2df(dat,s) |
44 |
} |
|
45 | ||
46 |
## Processing for final run |
|
47 | ! |
if (s == dat$control$nsims) { |
48 | ||
49 |
# Remove functions from control list |
|
50 | ! |
ftodel <- grep(".FUN", names(out$control), value = TRUE) |
51 | ! |
out$control[ftodel] <- NULL |
52 | ||
53 |
# Set column names for varying list elements |
|
54 | ! |
for (i in as.vector(which(lapply(out$epi, class) == "data.frame"))) { |
55 | ! |
colnames(out$epi[[i]]) <- paste0("sim", 1:dat$control$nsims) |
56 |
} |
|
57 | ||
58 |
} |
|
59 | ||
60 | ! |
return(out) |
61 |
} |
1 |
#' Sample value from Weibull distribution |
|
2 |
#' |
|
3 |
#' Function to sample values from a Weibull distribution with parameters shape |
|
4 |
#' and rate |
|
5 |
#' |
|
6 |
#' @param vecTimeSinceExp ? |
|
7 |
#' @param scale scale parameter value |
|
8 |
#' @param shape shape parameter value |
|
9 |
#' |
|
10 |
#' @return Value of coefficient of variation for vector |
|
11 |
#' @export |
|
12 |
cum_discr_si <- function(vecTimeSinceExp, scale, shape) { |
|
13 | ! |
vlen <- length(vecTimeSinceExp) |
14 | ! |
if (vlen > 0) { |
15 | ! |
probVec <- numeric(vlen) |
16 | ! |
for (p in 1:vlen) { |
17 | ! |
probVec[p] <- pweibull(vecTimeSinceExp[p], shape=shape, scale=scale) |
18 |
} |
|
19 |
} else { |
|
20 | ! |
probVec <- 0 |
21 |
} |
|
22 | ! |
return(probVec) |
23 |
} |
1 |
# Needed for recovery.FUN -------------------------------------------------------------- |
|
2 |
update_status <- function(rate, rand, active, status, label, state, at, expTime = NULL, prog.dist.scale = NULL, prog.dist.shape = NULL){ |
|
3 |
|
|
4 | 14x |
smp_sz <- 0 |
5 | 14x |
at_idx <- NULL |
6 |
|
|
7 | 14x |
idsElig <- which(active == 1 & status %in% label) |
8 | 14x |
nElig <- length(idsElig) |
9 |
|
|
10 | 14x |
if (nElig > 0) { |
11 | ||
12 | 10x |
if (rand) { |
13 | 6x |
vecProg <- which(stats::rbinom(nElig, 1, rate) == 1) |
14 | 6x |
if (length(vecProg) > 0) { |
15 | ! |
at_idx <- idsElig[vecProg] |
16 | ! |
smp_sz <- length(at_idx) |
17 | ! |
status[at_idx] <- state |
18 |
} |
|
19 |
}else{ |
|
20 | 4x |
vecTimeSinceExp <- at - expTime[idsElig] |
21 | 4x |
vecTimeSinceExp[is.na(vecTimeSinceExp)] <- 0 |
22 | 4x |
gammaRatesElig <- stats::pweibull(vecTimeSinceExp, prog.dist.shape, scale=prog.dist.scale) |
23 | 4x |
smp_sz <- round(sum(gammaRatesElig, na.rm=TRUE)) |
24 | 4x |
smp_prob <- gammaRatesElig |
25 |
|
|
26 | 4x |
if(smp_sz > 0){ |
27 | ! |
at_idx <- EpiModel::ssample(idsElig, smp_sz, prob = smp_prob) |
28 | ! |
status[at_idx] <- state |
29 |
} |
|
30 |
|
|
31 |
} |
|
32 |
|
|
33 |
} |
|
34 |
|
|
35 | 14x |
list(smp_sz, at_idx, status) |
36 |
} |
|
37 | ||
38 |
# Needed for departures.FUN ------------------------------------------------------------ |
|
39 |
update_active <- function(rate, rand, active, status, label){ |
|
40 |
|
|
41 | 66x |
smp_sz <- 0 |
42 | 66x |
act_idx <- NULL |
43 |
|
|
44 | 66x |
idsElig <- which(active == 1 & status == label) |
45 | 66x |
nElig <- length(idsElig) |
46 |
|
|
47 | 66x |
if (nElig > 0) { |
48 |
|
|
49 | 31x |
gElig <- rep(1, nElig) |
50 |
|
|
51 | 31x |
if (rand) { |
52 | 31x |
vec <- which(stats::rbinom(nElig, 1, rate) == 1) |
53 |
|
|
54 | 31x |
if (length(vec) > 0) { |
55 | 1x |
act_idx <- idsElig[vec] |
56 | 1x |
smp_sz <- length(act_idx) |
57 |
} |
|
58 |
} else { |
|
59 | ! |
smp_sz <- min(round(sum(rate[gElig == 1])), sum(gElig == 1)) |
60 | ! |
act_idx <- EpiModel::ssample(idsElig[gElig == 1], smp_sz) |
61 |
} |
|
62 |
} |
|
63 |
|
|
64 | 66x |
return(list(smp_sz, act_idx)) |
65 |
} |
1 |
#' Initialize ICM |
|
2 |
#' |
|
3 |
#' Function to initialize the icm |
|
4 |
#' |
|
5 |
#' @param param ICM parameters. |
|
6 |
#' @param init Initial value parameters. |
|
7 |
#' @param control Control parameters |
|
8 |
#' @param seed random seed for checking consistency with other versions. |
|
9 |
#' |
|
10 |
#' @return Updated dat |
|
11 |
#' @export |
|
12 |
initialize.FUN <- function(param, init, control, seed = NULL) { |
|
13 | 10x |
if(!is.null(seed)) set.seed(seed) |
14 |
|
|
15 |
## Master List for Data ## |
|
16 | 13x |
dat <- list() |
17 | 13x |
dat$param <- param |
18 | 13x |
dat$init <- init |
19 | 13x |
dat$control <- control |
20 |
|
|
21 |
# Set attributes |
|
22 | 13x |
dat$attr <- list() |
23 | 13x |
numeric.init <- init[which(sapply(init, class) == "numeric")] |
24 | 13x |
n <- do.call("sum", numeric.init) |
25 | 13x |
dat$attr$active <- rep(1, n) |
26 | 13x |
dat$attr$group <- rep(1, n) |
27 |
|
|
28 |
# Initialize status and infection time |
|
29 | 13x |
dat <- init_status.icm(dat) |
30 |
|
|
31 |
# Summary out list |
|
32 | 13x |
dat <- get_prev.icm(dat, at = 1) |
33 |
|
|
34 | 13x |
return(dat) |
35 |
} |
|
36 | ||
37 | ||
38 |
#' Initialize status of ICM |
|
39 |
#' |
|
40 |
#' Function to get the status of the initialized icm |
|
41 |
#' |
|
42 |
#' @param dat Object containing all data |
|
43 |
#' |
|
44 |
#' @return Updated dat |
|
45 |
#' @importFrom EpiModel ssample |
|
46 |
#' @export |
|
47 |
init_status.icm <- function(dat) { |
|
48 |
|
|
49 |
# Variables --------------------------------------------------------------- |
|
50 |
type <- dat$control$type |
|
51 |
group <- dat$attr$group |
|
52 |
nGroups <- dat$param$groups |
|
53 |
|
|
54 |
nG <- sum(group == 1) |
|
55 |
|
|
56 |
e.num <- dat$init$e.num |
|
57 |
i.num <- dat$init$i.num |
|
58 |
q.num <- dat$init$q.num |
|
59 |
h.num <- dat$init$h.num |
|
60 |
r.num <- dat$init$r.num |
|
61 |
f.num <- dat$init$f.num |
|
62 |
|
|
63 |
# Status ------------------------------------------------------------------ |
|
64 |
status <- rep("s", nG) |
|
65 |
status[sample(which(group == 1), size = i.num)] <- "i" |
|
66 |
if (type %in% c("SIR", "SEIR", "SEIQHR", "SEIQHRF")) { |
|
67 |
status[sample(which(group == 1 & status == "s"), size = r.num)] <- "r" |
|
68 |
} |
|
69 |
if (type %in% c("SEIR", "SEIQHR", "SEIQHRF")) { |
|
70 |
status[sample(which(group == 1 & status == "s"), size = e.num)] <- "e" |
|
71 |
} |
|
72 |
if (type %in% c("SEIQHR", "SEIQHRF")) { |
|
73 |
status[sample(which(group == 1 & status == "s"), size = q.num)] <- "q" |
|
74 |
status[sample(which(group == 1 & status == "s"), size = h.num)] <- "h" |
|
75 |
} |
|
76 |
if (type %in% c("SEIQHRF")) { |
|
77 |
status[sample(which(group == 1 & status == "s"), size = f.num)] <- "f" |
|
78 |
} |
|
79 |
|
|
80 |
dat$attr$status <- status |
|
81 |
n <- length(status) |
|
82 |
|
|
83 |
# leave exposure time uninitialised for now, and |
|
84 |
# just set to NA at start. |
|
85 |
|
|
86 |
# Exposure Time ---------------------------------------------------------- |
|
87 |
dat$attr$expTime <- rep(NA, n) |
|
88 |
|
|
89 |
# Infection Time ---------------------------------------------------------- |
|
90 |
infTime <- rep(NA, n) |
|
91 |
infTime[status == "i"] <- 1 |
|
92 |
dat$attr$infTime <- infTime |
|
93 |
|
|
94 |
# Recovery Time ---------------------------------------------------------- |
|
95 |
dat$attr$recovTime <- rep(NA, n) |
|
96 |
|
|
97 |
# Need for Hospitalisation Time ---------------------------------------------------------- |
|
98 |
dat$attr$hospTime <- rep(NA, n) |
|
99 |
|
|
100 |
# Quarantine Time ---------------------------------------------------------- |
|
101 |
dat$attr$quarTime <- rep(NA, n) |
|
102 |
|
|
103 |
# Hospital-need cessation Time ---------------------------------------------------------- |
|
104 |
dat$attr$dischTime <- rep(NA, n) |
|
105 |
|
|
106 |
# Case-fatality Time ---------------------------------------------------------- |
|
107 |
dat$attr$fatTime <- rep(NA, n) |
|
108 |
|
|
109 |
return(dat) |
|
110 |
} |
1 |
#' Generating time-dependent parameters |
|
2 |
#' |
|
3 |
#' Function to generate parameter values for the time range of the simulation |
|
4 |
#' that can change over time. |
|
5 |
#' |
|
6 |
#' @param nstep Number of time steps to generate parameter values for. |
|
7 |
#' @param vals List of parameter values to include over the nsteps. |
|
8 |
#' @param timing List of the step numbers at which to start changes, with the |
|
9 |
#' last number reflecting when to hit the last parameter value in vals. |
|
10 |
#' |
|
11 |
#' @return list of parameter values for length t |
|
12 |
#' |
|
13 |
#' @importFrom utils tail |
|
14 |
#' @export |
|
15 |
vary_param <- function(nstep, vals, timing) { |
|
16 | ||
17 | 1x |
stopifnot(length(vals) == length(timing)) |
18 | 1x |
y <- list() |
19 | ||
20 | 1x |
for(t in seq(1:nstep)){ |
21 | 90x |
if(t <= timing[1]){ # If before first jump set to val[1] |
22 | 7x |
y.t <- vals[1] |
23 | 83x |
}else if(t > utils::tail(timing, n=1)){ # If after last jump set to val[-1] |
24 | 76x |
y.t <- utils::tail(vals, n=1) |
25 |
}else{ # If intermediate step, calculate... |
|
26 | 7x |
for(j in (2:length(timing))){ |
27 | 7x |
if(t > timing[j - 1] & t <= timing[j]){ |
28 | 7x |
start <- vals[j - 1] |
29 | 7x |
end <- vals[j] |
30 | 7x |
start_t <- timing[j - 1] |
31 | 7x |
end_t <- timing[j] |
32 | 7x |
y.t <- start - (t-start_t)*(start - end)/(end_t - start_t) |
33 |
} |
|
34 |
} |
|
35 |
} |
|
36 | 90x |
y <- append(y, y.t) |
37 |
} |
|
38 | ||
39 | 1x |
return(unlist(y)) |
40 | ||
41 |
} |
1 |
#' summary function |
|
2 |
#' |
|
3 |
#' Function to extract mean and standard deviation of all compartments |
|
4 |
#' for SEIQHRF model |
|
5 |
#' |
|
6 |
#' @param object seiqhrf object, returned by \code{\link{seiqhrf}}. |
|
7 |
#' @param ... Additional parameters. |
|
8 |
#' |
|
9 |
#' @return A list of compartments, each compartment is a list of three: |
|
10 |
#' \itemize{ |
|
11 |
#' \item \code{mean}: (t-dimensional vector) Mean of simulated counts of the compartment; |
|
12 |
#' \item \code{sd}: (t-dimensional vector) Standard deviation of simulated counts of the compartment; |
|
13 |
#' \item \code{CI}: (t by 2 matrix) 95\% confidence intervals of simulated counts of the compartment (assumed Gaussian), |
|
14 |
#' \item \code{qntCI}: (t by 2 matrix) 95\% quantile confidence intervals, |
|
15 |
#' } |
|
16 |
#' where t is the total number of time points in the simulation (i.e. object$nstep). |
|
17 |
#' |
|
18 |
#' @importFrom stats setNames |
|
19 |
#' @export |
|
20 |
summary.seiqhrf <- function(object, ...){ |
|
21 |
|
|
22 | 2x |
nsteps <- object$control$nsteps |
23 | 2x |
c_names <- c("s.num", "e.num", "i.num", "q.num", "h.num", "r.num", "f.num") |
24 |
|
|
25 |
#### extract mean and sd |
|
26 | 2x |
df_mean <- as.data.frame(object, out = "mean") |
27 | 2x |
sub_cols <- names(df_mean)[ order(match(names(df_mean), c_names))][1:length(c_names)] |
28 | 2x |
df_mean <- df_mean[, sub_cols] |
29 |
|
|
30 | 2x |
df_sd <- as.data.frame(object, out = "sd")[, sub_cols] |
31 | ||
32 | 2x |
df_qnt_left <- as.data.frame(object, out = "qnt", qval = .025)[, sub_cols] |
33 | 2x |
df_qnt_right <- as.data.frame(object, out = "qnt", qval = .975)[, sub_cols] |
34 |
|
|
35 |
#### prepare result |
|
36 | 2x |
res <- stats::setNames(vector("list", length(c_names)), c_names) |
37 | ||
38 | 2x |
for(prev_no in c_names){ |
39 | ||
40 | 14x |
prev_mean <- df_mean[[prev_no]] |
41 | 14x |
prev_sd <- df_sd[[prev_no]] |
42 | 14x |
prev_ci <- cbind(prev_mean - 1.96*prev_sd, prev_mean + 1.96*prev_sd) |
43 | 14x |
prev_qci <- cbind(df_qnt_left[[prev_no]], df_qnt_right[[prev_no]]) |
44 | 14x |
names(prev_mean) <- c(1:nsteps) |
45 | 14x |
names(prev_sd) <- c(1:nsteps) |
46 | 14x |
rownames(prev_ci) <- c(1:nsteps) |
47 | 14x |
rownames(prev_qci) <- c(1:nsteps) |
48 |
|
|
49 | 14x |
res[[prev_no]]$"mean" <- prev_mean |
50 | 14x |
res[[prev_no]]$"sd" <- prev_sd |
51 | 14x |
res[[prev_no]]$"CI" <- prev_ci |
52 | 14x |
res[[prev_no]]$"qntCI" <- prev_qci |
53 |
|
|
54 |
} |
|
55 |
|
|
56 | 2x |
class(res) <- "summary.seiqhrf" |
57 | 2x |
res |
58 |
} |
|
59 |
1 |
#' Arrivals function |
|
2 |
#' |
|
3 |
#' Function to handle background demographics for the SEIQHRF model. |
|
4 |
#' Specifically arrivals (births and immigration). Uses the original EpiModel |
|
5 |
#' code currently. A replacement that implements modelling the arrival of |
|
6 |
#' infected individuals is under development -- but for now, all arrivals |
|
7 |
#' go into the S compartment.. |
|
8 |
#' |
|
9 |
#' @param dat Input data needed. |
|
10 |
#' @param at time point |
|
11 |
#' @param seed random seed for checking consistency |
|
12 |
#' |
|
13 |
#' @return Updated dat |
|
14 |
#' @export |
|
15 |
arrivals.FUN <- function(dat, at, seed = NULL) { |
|
16 |
|
|
17 | 10x |
if(!is.null(seed)) set.seed(seed) |
18 |
# Conditions -------------------------------------------------------------- |
|
19 | ! |
if (!dat$param$vital) return(dat) |
20 |
|
|
21 |
|
|
22 |
|
|
23 | 10x |
flare.inf.point <- dat$param$flare.inf.point |
24 | 10x |
flare.inf.num <- dat$param$flare.inf.num |
25 |
|
|
26 | 10x |
if(!(at %in% flare.inf.point)){ |
27 |
# Variables --------------------------------------------------------------- |
|
28 | 10x |
a.rate <- dat$param$a.rate |
29 | 10x |
a.prop.e <- dat$param$a.prop.e |
30 | 10x |
a.prop.i <- dat$param$a.prop.i |
31 | 10x |
a.prop.q <- dat$param$a.prop.q |
32 | 10x |
a.rand <- dat$control$a.rand |
33 | 10x |
nOld <- dat$epi$num[at - 1] |
34 |
# Process: partition arrivals into compartments ----------------------------------------------------------------- |
|
35 | 10x |
nArrivals <- ifelse(a.rand, sum(stats::rbinom(nOld, 1, a.rate)), round(nOld * a.rate)) |
36 | 10x |
nArrivals.e <- round(nArrivals * ifelse(length(a.prop.e) > 1, a.prop.e[at], a.prop.e)) |
37 | 10x |
nArrivals.i <- round(nArrivals * ifelse(length(a.prop.i) > 1, a.prop.i[at], a.prop.i)) |
38 | 10x |
nArrivals.q <- round(nArrivals * ifelse(length(a.prop.q) > 1, a.prop.q[at], a.prop.q)) |
39 | 10x |
totArrivals <- nArrivals.e + nArrivals.i + nArrivals.q |
40 | 10x |
dat$attr$status <- c(dat$attr$status, |
41 | 10x |
rep("e", nArrivals.e), |
42 | 10x |
rep("i", nArrivals.i), |
43 | 10x |
rep("q", nArrivals.q)) |
44 |
|
|
45 |
}else{ |
|
46 | ! |
nArrivals.i <- totArrivals <- flare.inf.num[which(flare.inf.point == at)] |
47 | ! |
dat$attr$status <- c(dat$attr$status, rep("i", nArrivals.i) ) |
48 | ! |
nArrivals.e <- nArrivals.q <- 0 |
49 |
} |
|
50 |
|
|
51 | 10x |
dat$attr$active <- c(dat$attr$active, rep(1, totArrivals)) |
52 | 10x |
dat$attr$group <- c(dat$attr$group, rep(1, totArrivals)) |
53 | 10x |
dat$attr$expTime <- c(dat$attr$expTime, rep(NA, totArrivals)) |
54 | 10x |
dat$attr$infTime <- c(dat$attr$infTime, rep(NA, totArrivals)) |
55 | 10x |
dat$attr$quarTime <- c(dat$attr$quarTime, rep(NA, totArrivals)) |
56 | 10x |
dat$attr$hospTime <- c(dat$attr$ihospTime, rep(NA, totArrivals)) |
57 | 10x |
dat$attr$recovTime <- c(dat$attr$recovTime, rep(NA, totArrivals)) |
58 | 10x |
dat$attr$fatTime <- c(dat$attr$fatTime, rep(NA, totArrivals)) |
59 |
|
|
60 |
|
|
61 |
# Output ------------------------------------------------------------------ |
|
62 | 10x |
if (at == 2) { |
63 | 10x |
dat$epi$a.flow <- c(0, totArrivals) |
64 | 10x |
dat$epi$a.e.flow <- c(0, nArrivals.e) |
65 | 10x |
dat$epi$a.i.flow <- c(0, nArrivals.i) |
66 | 10x |
dat$epi$a.q.flow <- c(0, nArrivals.q) |
67 |
} else { |
|
68 | ! |
dat$epi$a.flow[at] <- totArrivals |
69 | ! |
dat$epi$a.e.flow[at] <- c(0, nArrivals.e) |
70 | ! |
dat$epi$a.i.flow[at] <- c(0, nArrivals.i) |
71 | ! |
dat$epi$a.q.flow[at] <- c(0, nArrivals.q) |
72 |
} |
|
73 |
|
|
74 | 10x |
return(dat) |
75 |
} |
1 |
#' SEIQHRF Simulation Wrapper |
|
2 |
#' |
|
3 |
#' Wrapper to implement SEIQHRF model |
|
4 |
#' |
|
5 |
#' @param type Type of model: SI, SIR, SIS, SEIR, SEIQHR and |
|
6 |
#' SEIQHRF available, but only SEIQHRF is likely to work in the |
|
7 |
#' current version of the code. |
|
8 |
#' @param nsteps Number of time steps to solve the model over. This must be a |
|
9 |
#' positive integer. |
|
10 |
#' @param nsims Number of simulations to run. |
|
11 |
#' @param ncores Number of physical CPU cores used for parallel computation. |
|
12 |
#' @param prog.rand Method for progression from E compartment to I. If TRUE, |
|
13 |
#' random binomial draws at prog.rate, if FALSE, random draws from a |
|
14 |
#' Weibull distribution (yes, I know it should be a discrete Weibull |
|
15 |
#' distribution but it makes little difference and speed of computation |
|
16 |
#' matters), with parameters prog.dist.scale and prog.dist.shape |
|
17 |
#' @param rec.rand Method for recovery transition from I, Q or H to R. If TRUE, |
|
18 |
#' random binomial draws at rec.rate, if FALSE, random draws from a |
|
19 |
#' random draws from a Weibull distribution, with parameters |
|
20 |
#' rec.dist.scale and rec.dist.shape. |
|
21 |
#' @param arec.rand Method for recovery transition from E to R. If TRUE, |
|
22 |
#' random binomial draws at arec.rate, if FALSE, random draws from a |
|
23 |
#' random draws from a Weibull distribution, with parameters |
|
24 |
#' arec.dist.scale and arec.dist.shape. |
|
25 |
#' @param fat.rand Method for case fatality transition from H to F. If TRUE, |
|
26 |
#' random binomial draws at fat.rate.base, if FALSE, random sample with |
|
27 |
#' a sample fraction also given by fat.rate.base. However, if the |
|
28 |
#' current number of patients in the H (needs hospitalisation) |
|
29 |
#' compartment is above a hospital capacity level specified by |
|
30 |
#' hosp.cap, then the fatality rate is the mean of the base fatality |
|
31 |
#' rate weighted by the hospital capacity, plus a higher rate, |
|
32 |
#' specified by fat.rate.overcap, weighted by the balance of those |
|
33 |
#' requiring hospitalisation (but who can't be accommodated). By |
|
34 |
#' setting fat.rate.overcap higher, the effect of exceeding the |
|
35 |
#' capacity of the health care system can be simulated. There is also |
|
36 |
#' a coefficient fat.tcoeff for the fatality rates that increases them |
|
37 |
#' as a linear function of the number of days the individual has been |
|
38 |
#' in the H compartment. Use of the co-efficient better approximates |
|
39 |
#' the trapezoid survival time distribution typical of ICU patients. |
|
40 |
#' @param quar.rand Method for self-isolation transition from I to Q. If TRUE, |
|
41 |
#' random binomial draws at quar.rate, if FALSE, random sample with a |
|
42 |
#' sample fraction also given by `quar.rate. |
|
43 |
#' @param hosp.rand Method for transition from I or Q to H -- that is, from |
|
44 |
#' infectious or from self-isolated to requiring hospitalisation. If |
|
45 |
#' TRUE, random binomial draws at hosp.rate, if FALSE, random sample |
|
46 |
#' with a sample fraction also given by `hosp.rate. |
|
47 |
#' @param disch.rand Method for transition from H to R -- that is, from |
|
48 |
#' requiring hospitalisation to recovered. If TRUE, random binomial |
|
49 |
#' draws at disch.rate, if FALSE, random sample with a sample fraction |
|
50 |
#' also given by disch.rate. Note that the only way out of the H |
|
51 |
#' compartment is recovery or death. |
|
52 |
#' @param infection.FUN No, being infected with SARS-CoV2 is not fun. Rather |
|
53 |
#' this is the the name of the function to implement infection |
|
54 |
#' processes. Use the default. |
|
55 |
#' @param recovery.FUN Function to implement recovery processes. Use the |
|
56 |
#' default. |
|
57 |
#' @param departures.FUN Handles background demographics, specifically |
|
58 |
#' departures (deaths not due to the virus, and emigration). Use the |
|
59 |
#' default. |
|
60 |
#' @param arrivals.FUN Handles background demographics, specifically arrivals |
|
61 |
#' (births and immigration). Uses the original EpiModel code |
|
62 |
#' currently. A replacement that implements modelling the arrival of |
|
63 |
#' infected individuals is under development -- but for now, all |
|
64 |
#' arrivals go into the S compartment. |
|
65 |
#' @param get_prev.FUN Utility function that collects prevalence and transition |
|
66 |
#' time data from each run and stores it away in the simulation result |
|
67 |
#' object. Use the default. |
|
68 |
#' @param s.num Initial number of *S compartment individuals in |
|
69 |
#' the simulated population. An overall population of 10,000 is a good |
|
70 |
#' compromise. A set of models will still take several minutes or more |
|
71 |
#' to run, in parallel. |
|
72 |
#' @param e.num Initial number of E compartment individuals in |
|
73 |
#' the simulated population. |
|
74 |
#' @param i.num Initial number of I compartment individuals in |
|
75 |
#' the simulated population. |
|
76 |
#' @param q.num Initial number of Q compartment individuals in |
|
77 |
#' the simulated population. |
|
78 |
#' @param h.num Initial number of H compartment individuals in |
|
79 |
#' the simulated population. |
|
80 |
#' @param r.num Initial number of R compartment individuals in |
|
81 |
#' the simulated population. |
|
82 |
#' @param f.num Initial number of F compartment individuals in |
|
83 |
#' the simulated population. |
|
84 |
#' @param inf.prob.e Probability of passing on infection at each |
|
85 |
#' exposure event for interactions between infectious people in the E |
|
86 |
#' compartment and susceptible in S. Note the default is lower than for |
|
87 |
#' inf.prob.i reflecting the reduced infectivity of infected but |
|
88 |
#' asymptomatic people (the E compartment). Otherwise as for inf.exp.i. |
|
89 |
#' @param act.rate.e The number of exposure events (acts) between |
|
90 |
#' infectious individuals in the E compartment and susceptible individuals |
|
91 |
#' in the S compartment, per day. Otherwise as for act.rate.i. |
|
92 |
#' @param inf.prob.i Probability of passing on infection at each |
|
93 |
#' exposure event for interactions between infectious people in the I |
|
94 |
#' compartment and susceptible in S. Reducing inf.prob.i is equivalent to |
|
95 |
#' increasing hygiene measures, such as not putting hands in eyes, nose or |
|
96 |
#' moth, use of hand sanitisers, wearing masks by the infected, and so on. |
|
97 |
#' @param act.rate.i The number of exposure events (acts) between |
|
98 |
#' infectious individuals in the I compartment and susceptible individuals |
|
99 |
#' in the S compartment, per day. It's stochastic, so the rate is an |
|
100 |
#' average, some individuals may have more or less. Note that not every |
|
101 |
#' exposure event results in infection - that is governed by the inf.prob.i |
|
102 |
#' parameters (see below). Reducing act.rate.i is equivalent to increasing |
|
103 |
#' social distancing by people in the I compartment. |
|
104 |
#' @param inf.prob.q Probability of passing on infection at each |
|
105 |
#' exposure event for interactions between infectious people in the Q |
|
106 |
#' compartment and susceptible in S. Note the default is lower than for |
|
107 |
#' inf.prob.i reflecting the greater care that self-isolated individuals |
|
108 |
#' will, on average, take regarding hygiene measures, such as wearing masks, |
|
109 |
#' to limit spread to others. Otherwise as for inf.exp.i. |
|
110 |
#' @param act.rate.q The number of exposure events (acts) between |
|
111 |
#' infectious individuals in the Q compartment (isolated, self or otherwise) |
|
112 |
#' and susceptible individuals in the S compartment, per day. Note the much |
|
113 |
#' lower rate than for the I and E compartments, reflecting the much |
|
114 |
#' greater degree of social isolation for someone in (self-)isolation. The |
|
115 |
#' exposure event rate is not zero for this group, just much less. |
|
116 |
#' Otherwise as for act.rate.i. |
|
117 |
#' @param quar.rate Rate per day at which symptomatic (or tested |
|
118 |
#' positive), infected I compartment people enter self-isolation (Q |
|
119 |
#' compartment). Asymptomatic E compartment people can't enter |
|
120 |
#' self-isolation because they don't yet know they are infected. Default is |
|
121 |
#' a low rate reflecting low community awareness or compliance with |
|
122 |
#' self-isolation requirements or practices, but this can be tweaked when |
|
123 |
#' exploring scenarios. |
|
124 |
#' @param quar.dist.scale Scale parameter for Weibull distribution for |
|
125 |
#' recovery, see quar.rand for details. |
|
126 |
#' @param quar.dist.shape Shape parameter for Weibull distribution for |
|
127 |
#' recovery, see quar.rand for details. Read up on the Weibull distribution |
|
128 |
#' before changing the default. |
|
129 |
#' @param hosp.rate Rate per day at which symptomatic (or tested |
|
130 |
#' positive), infected I compartment people or self-isolated Q compartment |
|
131 |
#' people enter the state of requiring hospital care -- that is, become |
|
132 |
#' serious cases. A default rate of 1% per day with an average illness |
|
133 |
#' duration of about 10 days means a bit less than 10% of cases will |
|
134 |
#' require hospitalisation, which seems about right (but can be tweaked, |
|
135 |
#' of course). |
|
136 |
#' @param hosp.dist.scale Scale parameter for Weibull distribution for |
|
137 |
#' recovery, see hosp.rand for details. |
|
138 |
#' @param hosp.dist.shape Shape parameter for Weibull distribution for |
|
139 |
#' recovery, see hosp.rand for details. Read up on the Weibull distribution |
|
140 |
#' before changing the default. |
|
141 |
#' @param disch.rate Rate per day at which people needing |
|
142 |
#' hospitalisation recover. |
|
143 |
#' @param disch.dist.scale Scale parameter for Weibull distribution for |
|
144 |
#' recovery, see disch.rand for details. |
|
145 |
#' @param disch.dist.shape Shape parameter for Weibull distribution for |
|
146 |
#' recovery, see disch.rand for details. Read up on the Weibull distribution |
|
147 |
#' before changing the default. |
|
148 |
#' @param prog.rate Rate per day at which people who are infected |
|
149 |
#' but asymptomatic (E compartment) progress to becoming symptomatic (or |
|
150 |
#' test-positive), the I compartment. See prog.rand above for more details. |
|
151 |
#' @param prog.dist.scale Scale parameter for Weibull distribution |
|
152 |
#' for progression, see prog.rand for details. |
|
153 |
#' @param prog.dist.shape Shape parameter for Weibull distribution |
|
154 |
#' for progression, see prog.rand for details. Read up on the Weibull |
|
155 |
#' distribution before changing the default. |
|
156 |
#' @param rec.rate Rate per day at which people who are infected and |
|
157 |
#' symptomatic (I compartment) recover, thus entering the R compartment. |
|
158 |
#' See rec.rand above for more details. |
|
159 |
#' @param rec.dist.scale Scale parameter for Weibull distribution for |
|
160 |
#' recovery, see rec.rand for details. |
|
161 |
#' @param rec.dist.shape Shape parameter for Weibull distribution for |
|
162 |
#' recovery, see rec.rand for details. Read up on the Weibull distribution |
|
163 |
#' before changing the default. |
|
164 |
#' @param arec.rate Rate per day at which people who are exposed but asymptotic |
|
165 |
#' (E compartment) recover, thus entering the R compartment. |
|
166 |
#' See arec.rand above for more details. |
|
167 |
#' @param arec.dist.scale Scale parameter for Weibull distribution for |
|
168 |
#' recovery, see arec.rand for details. |
|
169 |
#' @param arec.dist.shape Shape parameter for Weibull distribution for |
|
170 |
#' recovery, see arec.rand for details. |
|
171 |
#' @param fat.rate.base Baseline mortality rate per day for people |
|
172 |
#' needing hospitalisation (deaths due to the virus). See fat.rand for more |
|
173 |
#' details. |
|
174 |
#' @param hosp.cap Number of available hospital beds for the modelled |
|
175 |
#' population. See fat.rand for more details. |
|
176 |
#' @param fat.rate.overcap Mortality rate per day for people needing |
|
177 |
#' hospitalisation but who can't get into hospital due to the hospitals |
|
178 |
#' being full (see hosp.cap and fat.rand). The default rate is twice that |
|
179 |
#' for those who do get into hospital. |
|
180 |
#' @param fat.tcoeff Time co-efficient for increasing mortality rate |
|
181 |
#' as time in the H compartment increases for each individual in it. See |
|
182 |
#' fat.rand for details. |
|
183 |
#' @param vital Enables demographics, that is, arrivals and |
|
184 |
#' departures, to and from the simulated population. |
|
185 |
#' @param a.rate Background demographic arrival rate. Currently all |
|
186 |
#' arrivals go into the S compartment, the default is approximately the |
|
187 |
#' daily birth rate for Australia. Will be extended to cover immigration in |
|
188 |
#' future versions. |
|
189 |
#' @param a.prop.e Arrivals proportion that goes to E (immigration). |
|
190 |
#' @param a.prop.i Arrivals proportion that goes to I (immigration). |
|
191 |
#' @param a.prop.q Arrivals proportion that goes to Q (immigration). |
|
192 |
#' @param ds.rate Background demographic departure (death not due to |
|
193 |
#' virus) rates. Defaults based on Australian crude death rates. Can be |
|
194 |
#' used to model emigration as well as deaths. |
|
195 |
#' @param de.rate Background demographic departure (death not due to |
|
196 |
#' virus) rates. Defaults based on Australian crude death rates. Can be |
|
197 |
#' used to model emigration as well as deaths. |
|
198 |
#' @param di.rate Background demographic departure (death not due to |
|
199 |
#' virus) rates. Defaults based on Australian crude death rates. Can be used |
|
200 |
#' to model emigration as well as deaths. |
|
201 |
#' @param dq.rate Background demographic departure (death not due to |
|
202 |
#' virus) rates. Defaults based on Australian crude death rates. Can be used |
|
203 |
#' to model emigration as well as deaths. |
|
204 |
#' @param dh.rate Background demographic departure (death not due to |
|
205 |
#' virus) rates. Defaults based on Australian crude death rates. Can be used |
|
206 |
#' to model emigration as well as deaths. |
|
207 |
#' @param dr.rate Background demographic departure (death not due to |
|
208 |
#' virus) rates. Defaults based on Australian crude death rates. Can be used |
|
209 |
#' to model emigration as well as deaths. |
|
210 |
#' @param out Summary function for the simulation runs. median is |
|
211 |
#' also available, or percentiles, see the EpiModel documentation. |
|
212 |
#' @param seed Random see for checking consistency. |
|
213 |
#' |
|
214 |
#' @return list with simulation and simulation dataframe |
|
215 |
#' @export |
|
216 |
simulate_seiqhrf <- function(type = "SEIQHRF", |
|
217 |
nsteps = 366, |
|
218 |
nsims = 8, |
|
219 |
ncores = 4, |
|
220 |
prog.rand = FALSE, |
|
221 |
quar.rand = TRUE, |
|
222 |
hosp.rand = TRUE, |
|
223 |
disch.rand = TRUE, |
|
224 |
rec.rand = FALSE, |
|
225 |
arec.rand = TRUE, |
|
226 |
fat.rand = TRUE, |
|
227 |
infection.FUN = 'infection.FUN', # infection.seiqhrf.icm, |
|
228 |
recovery.FUN = 'recovery.FUN', # progress.seiqhrf.icm, |
|
229 |
departures.FUN = 'departures.FUN', # departures.seiqhrf.icm, |
|
230 |
arrivals.FUN = 'arrivals.FUN', # arrivals.seiqhrf.icm, |
|
231 |
get_prev.FUN = 'get_prev.FUN', # get_prev.seiqhrf.icm, |
|
232 |
# init.icm params |
|
233 |
s.num = 9997, |
|
234 |
e.num=0, |
|
235 |
i.num = 3, |
|
236 |
q.num=0, |
|
237 |
h.num=0, |
|
238 |
r.num = 0, |
|
239 |
f.num = 0, |
|
240 |
# param.icm params |
|
241 |
inf.prob.e = 0.02, |
|
242 |
act.rate.e = 10, |
|
243 |
inf.prob.i = 0.05, |
|
244 |
act.rate.i = 10, |
|
245 |
inf.prob.q = 0.02, |
|
246 |
act.rate.q = 2.5, |
|
247 |
prog.rate = 1/10, |
|
248 |
quar.rate = 1/30, |
|
249 |
hosp.rate = 1/100, |
|
250 |
disch.rate = 1/15, |
|
251 |
rec.rate = 0.071, |
|
252 |
arec.rate = 0.05, |
|
253 |
prog.dist.scale = 5, |
|
254 |
prog.dist.shape = 1.5, |
|
255 |
quar.dist.scale = 1, |
|
256 |
quar.dist.shape = 1, |
|
257 |
hosp.dist.scale = 1, |
|
258 |
hosp.dist.shape = 1, |
|
259 |
disch.dist.scale = 1, |
|
260 |
disch.dist.shape = 1, |
|
261 |
rec.dist.scale = 35, |
|
262 |
rec.dist.shape = 1.5, |
|
263 |
arec.dist.scale = 35, |
|
264 |
arec.dist.shape = 1.5, |
|
265 |
fat.rate.base = 1/50, |
|
266 |
hosp.cap = 40, |
|
267 |
fat.rate.overcap = 1/25, |
|
268 |
fat.tcoeff = 0.5, |
|
269 |
vital = TRUE, |
|
270 |
a.rate = (10.5/365)/1000, |
|
271 |
a.prop.e = 0.01, |
|
272 |
a.prop.i = 0.001, |
|
273 |
a.prop.q = 0.01, |
|
274 |
ds.rate = (7/365)/1000, |
|
275 |
de.rate = (7/365)/1000, |
|
276 |
di.rate = (7/365)/1000, |
|
277 |
dq.rate = (7/365)/1000, |
|
278 |
dh.rate = (20/365)/1000, |
|
279 |
dr.rate = (7/365)/1000, |
|
280 |
out="mean", seed = NULL |
|
281 |
) { |
|
282 | ||
283 | ! |
control <- control_seiqhrf(type = type, |
284 | ! |
nsteps = nsteps, |
285 | ! |
nsims = nsims, |
286 | ! |
ncores = ncores, |
287 | ! |
prog.rand = prog.rand, |
288 | ! |
quar.rand = quar.rand, |
289 | ! |
hosp.rand = hosp.rand, |
290 | ! |
disch.rand = disch.rand, |
291 | ! |
rec.rand = rec.rand, |
292 | ! |
arec.rand = arec.rand, |
293 | ! |
infection.FUN = infection.FUN, |
294 | ! |
recovery.FUN = recovery.FUN, |
295 | ! |
arrivals.FUN = arrivals.FUN, |
296 | ! |
departures.FUN = departures.FUN, |
297 | ! |
get_prev.FUN = get_prev.FUN) |
298 | ||
299 | ! |
init <- init_seiqhrf(s.num = s.num, |
300 | ! |
e.num = e.num, |
301 | ! |
i.num = i.num, |
302 | ! |
q.num = q.num, |
303 | ! |
h.num = h.num, |
304 | ! |
r.num = r.num, |
305 | ! |
f.num = f.num) |
306 | ||
307 | ! |
param <- param_seiqhrf(inf.prob.e = inf.prob.e, |
308 | ! |
act.rate.e = act.rate.e, |
309 | ! |
inf.prob.i = inf.prob.i, |
310 | ! |
act.rate.i = act.rate.i, |
311 | ! |
inf.prob.q = inf.prob.q, |
312 | ! |
act.rate.q = act.rate.q, |
313 | ! |
prog.rate = prog.rate, |
314 | ! |
quar.rate = quar.rate, |
315 | ! |
hosp.rate = hosp.rate, |
316 | ! |
disch.rate = disch.rate, |
317 | ! |
rec.rate = rec.rate, |
318 | ! |
arec.rate = arec.rate, |
319 | ! |
prog.dist.scale = prog.dist.scale, |
320 | ! |
prog.dist.shape = prog.dist.shape, |
321 | ! |
quar.dist.scale = quar.dist.scale, |
322 | ! |
quar.dist.shape = quar.dist.shape, |
323 | ! |
hosp.dist.scale = hosp.dist.scale, |
324 | ! |
hosp.dist.shape = hosp.dist.shape, |
325 | ! |
disch.dist.scale = disch.dist.scale, |
326 | ! |
disch.dist.shape = disch.dist.shape, |
327 | ! |
rec.dist.scale = rec.dist.scale, |
328 | ! |
rec.dist.shape = rec.dist.shape, |
329 | ! |
arec.dist.scale = arec.dist.scale, |
330 | ! |
arec.dist.shape = arec.dist.shape, |
331 | ! |
fat.rate.base = fat.rate.base, |
332 | ! |
hosp.cap = hosp.cap, |
333 | ! |
fat.rate.overcap = fat.rate.overcap, |
334 | ! |
fat.tcoeff = fat.tcoeff, |
335 | ! |
vital = vital, |
336 | ! |
a.rate = a.rate, |
337 | ! |
a.prop.e = a.prop.e, |
338 | ! |
a.prop.i = a.prop.i, |
339 | ! |
a.prop.q = a.prop.q, |
340 | ! |
ds.rate = ds.rate, |
341 | ! |
de.rate = de.rate, |
342 | ! |
di.rate = di.rate, |
343 | ! |
dq.rate = dq.rate, |
344 | ! |
dh.rate = dh.rate, |
345 | ! |
dr.rate = dr.rate) |
346 | ||
347 | ! |
sim <- seiqhrf(param = param, init = init, control = control, seed = seed) |
348 | ! |
sim_df <- as.data.frame(sim, out=out) |
349 | ! |
class(sim) <- "icm" |
350 | ||
351 | ! |
return(list(sim=sim, df=sim_df)) |
352 |
} |