forked from epiforecasts/ringbp
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathoutbreak_step.R
169 lines (146 loc) · 6.95 KB
/
outbreak_step.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
#' Move forward one generation in the branching process
#'
#' @author Joel Hellewell
#'
#' @param case_data a `data.table`: cases in outbreak so far; initially
#' generated by [outbreak_setup()]
#' @inheritParams outbreak_model
#' @param incfn a `function` that samples from incubation period Weibull
#' distribution; generated using [dist_setup()]
#' @param delayfn a `function` that samples from the onset-to-hospitalisation
#' delay Weibull distribution; generated using [dist_setup()]
#'
#' @importFrom data.table data.table rbindlist
#' @importFrom purrr map2 map2_dbl map_lgl
#' @importFrom stats rbinom
#'
#' @return A `list` with three elements:
#' 1. `$cases`: a `data.table` with case data
#' 2. `$effective_r0`: a `numeric` with the effective reproduction number
#' 3. `$cases_in_gen`: a `numeric` with the number of new cases in that
#' generation
#' @export
#'
#' @examples
#'
#'\dontrun{
#' # incubation period sampling function
#' incfn <- dist_setup(dist_shape = 2.322737,dist_scale = 6.492272)
#' # delay distribution sampling function
#' delayfn <- dist_setup(1.651524, 4.287786)
#' # generate initial cases
#' case_data <- outbreak_setup(num.initial.cases = 5,incfn,delayfn,k=1.95,prop.asym=0)
#' # generate next generation of cases
#' case_data <- outbreak_step(case_data = case_data,
#' disp.iso = 1,
#' disp.com = 0.16,
#' r0isolated = 0,
#' r0subclin = 1.25,
#' disp.subclin = 0.16,
#' r0community = 2.5,
#' prop.asym = 0,
#' incfn = incfn,
#' delayfn = delayfn,
#' prop.ascertain = 0,
#' k = 1.95,
#' quarantine = FALSE)[[1]]
#'}
outbreak_step <- function(case_data = NULL, disp.iso = NULL, disp.com = NULL,
r0isolated = NULL, r0community = NULL,
prop.asym = NULL, incfn = NULL, delayfn = NULL,
prop.ascertain = NULL, k = NULL, quarantine = NULL,
r0subclin = NULL, disp.subclin = NULL) {
# For each case in case_data, draw new_cases from a negative binomial distribution
# with an R0 and dispersion dependent on if isolated=TRUE
case_data[, new_cases := purrr::map2_dbl(
ifelse(vect_isTRUE(isolated),
disp.iso,
ifelse(vect_isTRUE(asym), disp.subclin, disp.com)),
ifelse(vect_isTRUE(isolated),
r0isolated,
ifelse(vect_isTRUE(asym),r0subclin, r0community)),
~ rnbinom(1, size = .x, mu = .y))
]
# Select cases that have generated any new cases
new_case_data <- case_data[new_cases > 0]
# The total new cases generated
total_new_cases <- case_data[, sum(new_cases), ]
# If no new cases drawn, outbreak is over so return case_data
if (total_new_cases == 0) {
# If everyone is isolated it means that either control has worked or everyone has had a chance to infect but didn't
case_data$isolated <- TRUE
effective_r0 <- 0
cases_in_gen <- 0
out <- list(case_data, effective_r0, cases_in_gen)
names(out) <- c("cases", "effective_r0", "cases_in_gen")
return(out)
}
# Compile a data.table for all new cases, new_cases is the amount of people that each infector has infected
inc_samples <- incfn(total_new_cases)
prob_samples <- data.table(
# time when new cases were exposed, a draw from serial interval based on infector's onset
exposure = unlist(purrr::map2(new_case_data$new_cases, new_case_data$onset,
function(x, y) {
inf_fn(rep(y, x), k)
})),
# records the infector of each new person
infector = unlist(purrr::map2(new_case_data$caseid, new_case_data$new_cases,
function(x, y) {
rep(as.integer(x), as.integer(y))
})),
# records when infector was isolated
infector_iso_time = unlist(purrr::map2(new_case_data$isolated_time, new_case_data$new_cases,
function(x, y) {
rep(x, as.integer(y))
})),
# records if infector asymptomatic
infector_asym = unlist(purrr::map2(new_case_data$asym, new_case_data$new_cases,
function(x, y) {
rep(x, y)
})),
# draws a sample to see if this person is asymptomatic
asym = as.logical(rbinom(n = total_new_cases, 1, prob = prop.asym)),
# draws a sample to see if this person is traced
missed = as.logical(rbinom(n = total_new_cases, 1, prob = 1 - prop.ascertain)),
# sample from the incubation period for each new person
incubfn_sample = inc_samples,
isolated = FALSE,
new_cases = NA
)
prob_samples <- prob_samples[exposure < infector_iso_time][, # filter out new cases prevented by isolation
`:=`(# onset of new case is exposure + incubation period sample
onset = exposure + incubfn_sample)]
# cases whose parents are asymptomatic are automatically missed
prob_samples$missed[vect_isTRUE(prob_samples$infector_asym)] <- TRUE
# If you are asymptomatic, your isolation time is Inf
prob_samples[, isolated_time := ifelse(vect_isTRUE(asym), Inf,
# If you are not asymptomatic, but you are missed,
# you are isolated at your symptom onset
ifelse(vect_isTRUE(missed), onset + delayfn(1),
# If you are not asymptomatic and you are traced,
# you are isolated at max(onset,infector isolation time) # max(onset,infector_iso_time)
ifelse(!vect_isTRUE(rep(quarantine, total_new_cases)),
pmin(onset + delayfn(1), pmax(onset, infector_iso_time)),
infector_iso_time)))]
# Chop out unneeded sample columns
prob_samples[, c("incubfn_sample", "infector_iso_time", "infector_asym") := NULL]
# Set new case ids for new people
prob_samples$caseid <- (nrow(case_data) + 1):(nrow(case_data) + nrow(prob_samples))
## Number of new cases
cases_in_gen <- nrow(prob_samples)
## Estimate the effective r0
effective_r0 <- nrow(prob_samples) / nrow(case_data[!vect_isTRUE(case_data$isolated)])
# Everyone in case_data so far has had their chance to infect and are therefore considered isolated
case_data$isolated <- TRUE
# bind original cases + new secondary cases
case_data <- data.table::rbindlist(list(case_data, prob_samples),
use.names = TRUE)
# Return
out <- list(case_data, effective_r0, cases_in_gen)
names(out) <- c("cases", "effective_r0", "cases_in_gen")
return(out)
}
# A vectorised version of isTRUE
vect_isTRUE <- function(x) {
purrr::map_lgl(x, isTRUE)
}