-
Notifications
You must be signed in to change notification settings - Fork 0
/
est_relaxation_rt.R
88 lines (73 loc) · 2.38 KB
/
est_relaxation_rt.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
suppressPackageStartupMessages({
require(data.table)
require(EpiNow2)
require(qs)
})
#' TODO best way to get these arguments in? dates could come from intervention timing
.debug <- c("~/Dropbox/SA2UK", "ZAF", "2020-06-01", "2020-10-15")
.args <- if (interactive()) sprintf(c(
"%s/inputs/epi_data.rds",
"%s/inputs/pops/%s.rds",
"%s/inputs/covidm_fit_yu.qs",
.debug[3],
.debug[4],
.debug[2],
"%s/outputs/relaxation/%s.rds"
), .debug[1], .debug[2]) else commandArgs(trailingOnly = TRUE)
est.window <- 14
tariso <- tail(.args, 2)[1]
window.start <- as.Date(.args[2])
window.end <- as.Date(.args[3])
case.dt <- setkey(
readRDS(.args[1])[iso3==tariso, .(date, confirm=cases)], date
)[between(date, window.start-est.window, window.end+est.window*2)]
params <- readRDS(.args[4])
yu_fits <- qread(.args[5])[order(ll)]
yu_fits[, eqs := (1:.N)/.N ]
#' using the median yu fits
medyu <- yu_fits[which.max(eqs > 0.5)]
yref <- unname(as.matrix(medyu[, .SD, .SDcols = grep("y_",colnames(medyu))]))
uref <- unname(as.matrix(medyu[, .SD, .SDcols = grep("u_",colnames(medyu))]))
ys <- rep(yref[1, ], each = 2)
us <- rep(uref[1, ], each = 2)
params$pop <- lapply(
params$pop,
function(x){
x$y <- ys
x$u <- us
return(x)
}
)
load("NGM.rda")
# Set up example generation time
generation_time <- as.list(
EpiNow2::generation_times[disease == "SARS-CoV-2",
.(mean, mean_sd, sd, sd_sd, max=30)
])
tarmcv <- generation_time$mean_sd/generation_time$mean
tarscv <- generation_time$sd_sd/generation_time$sd
tarcv <- generation_time$sd/generation_time$mean
generation_time$mean <- unname(cm_generation_time(params))
generation_time$mean_sd <- generation_time$mean * tarmcv
generation_time$sd <- generation_time$mean * tarcv
generation_time$sd_sd <- generation_time$sd * tarscv
# Set delays between infection and case report
# (any number of delays can be specifed here)
incubation_period <- as.list(
EpiNow2::incubation_periods[disease == "SARS-CoV-2",
.(mean, mean_sd, sd, sd_sd, max=30)
])
re.est <- estimate_infections(
reported_cases = case.dt,
generation_time = generation_time,
delays = delay_opts(incubation_period),
stan = stan_opts(
samples = 4e3,
warmup = 200,
cores = 4,
control = list(adapt_delta = 0.9)
),
horizon = 0, CrIs = c(0.5, 0.95),
verbose = TRUE
)$summarised[variable == "R"][between(date, window.start, window.end)]
saveRDS(re.est, tail(.args, 1))