Skip to content

Commit ea0e1b0

Browse files
kkmannkkmann
authored
Merge pull request #2 from kkmann/extended-sensitivity-analyses
Extended sensitivity analyses
2 parents f27cd34 + 17f9340 commit ea0e1b0

20 files changed

+682
-185
lines changed

‎.gitignore

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -10,3 +10,5 @@ figures
1010
*_cache
1111
*_files
1212
*.html
13+
14+
.rcache

‎.julia/Manifest.toml

Lines changed: 10 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -110,9 +110,9 @@ version = "0.26.3"
110110

111111
[[FillArrays]]
112112
deps = ["LinearAlgebra", "Random", "SparseArrays"]
113-
git-tree-sha1 = "4705cc4e212c3c978c60b1b18118ec49b4d731fd"
113+
git-tree-sha1 = "dd4ab4257c257532003eb9836eea07473fcc732e"
114114
uuid = "1a297f60-69ca-5386-bcde-b61e274b549b"
115-
version = "0.11.5"
115+
version = "0.11.6"
116116

117117
[[Formatting]]
118118
deps = ["Printf"]
@@ -270,10 +270,10 @@ uuid = "79098fc4-a85e-5d69-aa6a-4863f24498fa"
270270
version = "0.6.1"
271271

272272
[[Rmath_jll]]
273-
deps = ["Libdl", "Pkg"]
274-
git-tree-sha1 = "d76185aa1f421306dec73c057aa384bad74188f0"
273+
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
274+
git-tree-sha1 = "1b7bf41258f6c5c9c31df8c1ba34c1fc88674957"
275275
uuid = "f50d1b31-88e8-58de-be2c-1cc44531875f"
276-
version = "0.2.2+1"
276+
version = "0.2.2+2"
277277

278278
[[SHA]]
279279
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"
@@ -344,9 +344,9 @@ version = "1.0.0"
344344

345345
[[Tables]]
346346
deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "LinearAlgebra", "TableTraits", "Test"]
347-
git-tree-sha1 = "a716dde43d57fa537a19058d044b495301ba6565"
347+
git-tree-sha1 = "f03fc113290ee7726b173fc7ea661260d204b3f2"
348348
uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
349-
version = "1.3.2"
349+
version = "1.4.0"
350350

351351
[[Test]]
352352
deps = ["Distributed", "InteractiveUtils", "Logging", "Random"]
@@ -366,9 +366,9 @@ uuid = "efce3f68-66dc-5838-9240-27a6d6f5f9b6"
366366
version = "0.5.3"
367367

368368
[[cov19sim]]
369-
deps = ["DataFrames", "DataFramesMeta", "Distributed", "Distributions", "Documenter", "Interpolations", "Random", "Statistics", "Test", "UUIDs"]
370-
git-tree-sha1 = "d200e3affaab6a4b371ef354dcb94797992c56c0"
371-
repo-rev = "main"
369+
deps = ["DataFrames", "DataFramesMeta", "Distributed", "Distributions", "Documenter", "Interpolations", "LinearAlgebra", "Random", "Statistics", "Test", "UUIDs"]
370+
git-tree-sha1 = "88272b4ac28a5da90e55f97df0d5ca143faa39a6"
371+
repo-rev = "re-introduce-random-effect-for-test-sensitivity"
372372
repo-url = "https://github.com/kkmann/cov19sim"
373373
uuid = "7070818b-bc4a-4efc-9908-1d2131e194b1"
374374
version = "0.1.0"

‎code/evaluate_performance.jl

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -8,7 +8,7 @@ function f(school, pr_external_infections, days; n = 1L)
88
end
99
step!(x)
1010
end
11-
return evaluate(x)
11+
return evaluate(x; tests = ["pcr", "lfd"])
1212
end
1313
return vcat(pmap(g, [resample(school) for i = 1:n])...)
1414
end

‎code/scenario.R

Lines changed: 127 additions & 52 deletions
Original file line numberDiff line numberDiff line change
@@ -13,12 +13,17 @@ scenario <- function(
1313
pcr_lod = 300.0,
1414
pcr_sens = .975,
1515
pcr_spec = 1.0,
16-
gamma = 0.0,
17-
eta = 1.0,
1816
lfd_spec = 0.998,
17+
lfd_ranef = 0.0,
1918
ar_window = 3L,
2019
ar_coefficient = 0.0,
2120
days = as.integer(6*7),
21+
scale = 0.0,
22+
l = 2.5,
23+
df = 3.0,
24+
gamma_min = 0.0, gamma_max = 0.1,
25+
rzero = 3.0,
26+
mean_sensitivity = 0.6,
2227
...
2328
) {
2429
res <- as.list(environment())
@@ -56,80 +61,143 @@ scenario <- function(
5661

5762
res
5863
}
59-
60-
sensitivity <- function(vl, params) {
61-
lfd <- julia_call("LogRegTest", "lfd", params$lfd_slope, params$lfd_intercept, params$lfd_spec, need_return = "Julia")
62-
julia_call("sensitivity.", lfd, vl, need_return = "R")
63-
}
6464

6565
n_class <- function(params) with(params, n_bubble * bubbles_per_class)
6666
n_school <- function(params) with(params, classes_per_school * n_class(params))
6767
n_weekly_infections <- function(params) with(params, 7*n_school(params)*pr_external_infection)
6868

6969
schooldays <- function(params) with(params, sum(((1:days) %% 7) %in% 1:5))
7070

71+
disease_model <- function(params, gamma = NULL) {
72+
if (is.null(gamma))
73+
gamma <- gamma(params$rzero, params)
74+
dm <- julia_call("LarremoreModel",
75+
gamma, frac_symptomatic = params$frac_symptomatic,
76+
l10vl_clearance = log10(params$lli),
77+
need_return = "Julia"
78+
)
79+
if (params$scale > 1) {
80+
dm <- julia_call("HeavyTailsModel",
81+
dm,
82+
l = params$l,
83+
scale = params$scale,
84+
df = params$df
85+
)
86+
}
87+
return(dm)
88+
}
89+
7190
pcr <- function(params) with(params, julia_call("FixedTest", "pcr", pcr_sens, pcr_spec, lod = pcr_lod, need_return = "Julia"))
7291

73-
lfd <- function(params) {
74-
with(params,
75-
julia_call("LogRegTest", "lfd",
76-
eta*params$lfd_slope, params$lfd_intercept, lfd_spec,
77-
ar_window = ar_window, ar_coefficient = ar_coefficient
78-
)
92+
lfd <- function(params, eta = NULL) {
93+
if (is.null(eta))
94+
eta <- eta(params$mean_sensitivity, params)
95+
julia_call("LogRegTest", "lfd",
96+
eta*params$lfd_slope, params$lfd_intercept, params$lfd_spec,
97+
ar_window = params$ar_window, ar_coefficient = params$ar_coefficient,
98+
ranef = params$lfd_ranef
7999
)
80100
}
81101

82-
fit_rzero_ <- function(params, gamma_min = 0, gamma_max = .1) {
102+
fit_rzero <- memoise::memoise(function(
103+
n_bubble, bubbles_per_class, classes_per_school, pr_meet_class, pr_meet_school,
104+
frac_symptomatic, pr_noncovid_symptoms, l, scale, df, lli, gamma_min, gamma_max
105+
) {
106+
get_school <- function(gamma) {
107+
julia_call("school",
108+
n_bubble, bubbles_per_class, classes_per_school, pr_meet_class, pr_meet_school,
109+
gamma,
110+
frac_symptomatic, pr_noncovid_symptoms, l, scale, df,
111+
Inf, 1, lli, # test specific values do not affect r0
112+
julia_call("DoNothing")
113+
)
114+
}
83115
tbl_rzero <- tibble(
84116
gamma = seq(gamma_min, gamma_max, length.out = 100),
85-
R = map(gamma, ~{
86-
params$gamma <- .
87-
sample_rzero(do.call(school, params), n = 10L)
88-
})
117+
R = map(gamma, ~sample_rzero(get_school(.), n = 10L))
89118
) %>%
90119
unnest(R)
91120
fit <- lm(formula = R ~ gamma - 1, data = tbl_rzero)
92121
list(fit = fit, data = tbl_rzero)
93-
}
94-
fit_rzero <- memoise::memoise(fit_rzero_)
95-
rzero <- function(gamma, params, gamma_min = 0, gamma_max = .1) {
96-
fit <- fit_rzero(params, gamma_min, gamma_max)$fit
122+
}, cache = mcache)
123+
rzero <- memoise::memoise(function(gamma, params) {
124+
fit <- with(params, fit_rzero(
125+
n_bubble, bubbles_per_class, classes_per_school, pr_meet_class, pr_meet_school,
126+
frac_symptomatic, pr_noncovid_symptoms, l, scale, df, lli, gamma_min, gamma_max
127+
))$fit
97128
as.numeric(predict(fit, newdata = tibble(gamma = gamma), type = "response"))
98-
}
99-
gamma <- function(R, params, gamma_min = 0, gamma_max = .1) {
100-
as.numeric(uniroot(function(x) rzero(x, params, gamma_min, gamma_max) - R, interval = c(gamma_min, gamma_max))$root)
101-
}
129+
}, cache = mcache)
130+
gamma <- memoise::memoise(function(R, params) {
131+
uniroot(
132+
function(x) rzero(x, params) - R,
133+
interval = c(params$gamma_min, params$gamma_max)
134+
)$root %>%
135+
as.numeric()
136+
}, cache = mcache)
102137

103-
sample_presymptomatic_vl <- function(params, n = 1e5) {
104-
dm <- julia_call("LarremoreModel", params$gamma, frac_symptomatic = params$frac_symptomatic, need_return = "Julia")
105-
individuals <- julia_call("Individual.", dm, rep(params$pr_noncovid_symptoms, n), need_return = "Julia")
138+
sample_presymptomatic_vl <- function(
139+
disease_model,
140+
pr_noncovid_symptoms,
141+
pcr_lod
142+
) {
143+
individuals <- julia_call("Individual.",
144+
disease_model,
145+
rep(pr_noncovid_symptoms, 1e4),
146+
need_return = "Julia"
147+
)
106148
julia_call("infect!.", individuals, need_return = "Julia")
107149
julia_call("steps!.", individuals, 21L, need_return = "Julia")
150+
tbl_u <- tibble(
151+
uuid = julia_call("string.",
152+
julia_call("getproperty.",
153+
individuals, julia_call("Symbol", "uuid")
154+
)
155+
),
156+
u = julia_call("getproperty.",
157+
individuals, julia_call("Symbol", "u_sensitivity", need_return = "Julia")
158+
)
159+
)
108160
julia_call("get_status_logs", individuals) %>%
109161
as_tibble() %>%
162+
left_join(tbl_u, by = "uuid") %>%
110163
arrange(uuid, day) %>%
111164
group_by(uuid) %>%
112165
filter(
113-
row_number() >= which(viral_load > params$pcr_lod)[1],
114-
row_number() < which(symptomatic)[1]
166+
row_number() >= which(viral_load > pcr_lod)[1],
167+
row_number() < which(symptomatic)[1] %>% {if_else(is.na(.), Inf, as.numeric(.))}
115168
) %>%
116169
sample_n(1) %>%
117170
ungroup() %>%
118-
select(uuid, day, viral_load)
171+
select(
172+
uuid,
173+
day,
174+
u,
175+
viral_load
176+
)
119177
}
120-
sample_presymptomatic_vl_ <- memoise::memoise(sample_presymptomatic_vl)
121-
mean_sensitivity <- function(params, eta, ...) {
122-
sample_presymptomatic_vl_(params, ...) %>%
123-
mutate(
124-
sensitivity = sensitivity(viral_load^eta, params)
178+
sample_presymptomatic_vl_ <- memoise::memoise(sample_presymptomatic_vl, cache = mcache)
179+
mean_sensitivity <- memoise::memoise(function(eta, params) {
180+
lfd <- lfd(params, eta = eta)
181+
sample_presymptomatic_vl_(
182+
disease_model(params, gamma = 0),
183+
params$pr_noncovid_symptoms,
184+
params$pcr_lod
125185
) %>%
126-
pull(sensitivity) %>%
127-
summary() %>%
186+
mutate(
187+
sensitivity = map2_dbl(viral_load, u,
188+
~julia_call("get_probability_positive", lfd, ..1, u = ..2)
189+
)
190+
) %>%
191+
pull(sensitivity) %>%
128192
mean()
129-
}
130-
eta <- function(params, target, ...) {
131-
as.numeric(uniroot(function(x) mean_sensitivity(params, x, ...) - target, interval = c(0, 5))$root)
132-
}
193+
}, cache = mcache)
194+
eta <- memoise::memoise(function(target, params) {
195+
uniroot(
196+
function(x) mean_sensitivity(x, params) - target,
197+
interval = c(0, 3)
198+
)$root %>%
199+
as.numeric()
200+
}, cache = mcache)
133201

134202
expand_scenario <- function(params = scenario(), ...) {
135203
tbl <- if (length(list(...)) == 0) {
@@ -163,20 +231,27 @@ evaluate_performance_ <- function(params, policy, n = 1L) {
163231
`% schooldays missed (cumulative)` = workdays_missed/n_school(params)/schooldays(params)
164232
)
165233
}
166-
evaluate_performance <- function(policies, params = scenario(), n = if (!is.null(n_resample)) {n_resample} else {25L}, ...) {
167-
gamma2rs <- memoise::memoise(function(gamma) round(rzero(gamma, params), 1) )
168-
eta2mean_sensitivity <- memoise::memoise(function(eta) round(mean_sensitivity(params, eta), 2) )
234+
evaluate_performance_mem <- memoise::memoise(evaluate_performance_, cache = mcache)
235+
evaluate_performance <- function(
236+
policies,
237+
params = scenario(),
238+
.gamma_min = 0.0, .gamma_max = 0.1,
239+
n = if (!is.null(n_resample)) {n_resample} else {25L},
240+
...
241+
) {
169242
expand_scenario(params, ...) %>%
170243
expand_grid(tibble(policy = policies)) %>%
171244
mutate(
172245
policy_name = names(policy),
173-
results = map2(params, policy, evaluate_performance_, n)
246+
results = map2(params, policy, evaluate_performance_mem, n)
174247
) %>%
175248
unnest(results) %>%
176-
select(-params) %>%
177249
mutate(
178-
policy_name = factor(policy_name, levels = names(lst_policies)),
179-
Rs = map_dbl(.data$gamma, gamma2rs),
180-
`mean sensitivity` = map_dbl(.data$eta, eta2mean_sensitivity)
181-
)
250+
policy_name = factor(policy_name, levels = names(lst_policies))
251+
) %>%
252+
rename(
253+
Rs = rzero,
254+
`mean sensitivity` = mean_sensitivity
255+
) %>%
256+
select(-params)
182257
}

‎code/school.R

Lines changed: 4 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,12 +1,14 @@
11
julia_eval('@everywhere include("code/school.jl")')
22

3-
school <- function(policy = julia_call("DoNothing", need_return = "Julia"), ...){
3+
school <- function(policy = julia_call("DoNothing", need_return = "Julia"), gamma = NULL, ...){
44
params <- scenario(...)
5+
if (is.null(gamma))
6+
gamma <- gamma(params$rzero, params)
57
with(params,
68
julia_call("school",
79
n_bubble, bubbles_per_class, classes_per_school, pr_meet_class, pr_meet_school,
810
gamma,
9-
frac_symptomatic, pr_noncovid_symptoms,
11+
frac_symptomatic, pr_noncovid_symptoms, l, scale, df,
1012
a, b,
1113
lli,
1214
policy

‎code/school.jl

Lines changed: 6 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,9 +1,14 @@
11
function school(
22
n_per_bubble, bubbles_per_class, m_classes, pr_meet_class, pr_meet_school,
3-
gamma, frac_symptomatic, pr_noncov_symptoms, a, b, lli,
3+
gamma,
4+
frac_symptomatic, pr_noncov_symptoms, l, scale, df,
5+
a, b, lli,
46
policy
57
)
68
dm = LarremoreModel(gamma; frac_symptomatic = frac_symptomatic, l10vl_clearance = log10(lli))
9+
if scale > 0
10+
dm = HeavyTailsModel(dm; l = l, scale = scale, df = df)
11+
end
712
ThreeLevelPopulation(
813
n_per_bubble = n_per_bubble,
914
bubbles_per_class = bubbles_per_class,

‎code/setup.jl

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -3,7 +3,7 @@ Pkg.activate(".julia")
33
Pkg.instantiate()
44
Pkg.precompile()
55

6-
using cov19sim, DataFrames, Distributions
6+
@time using cov19sim, DataFrames, Distributions, Random
77

88
import Random.seed!
99
seed!(42)

‎plots/autocorrelation.R

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
get_ar <- function(params, seed = 42L, nrsmpl = 1e4) {
22
julia_call("Random.seed!", as.integer(seed))
33
lfd <- lfd(params)
4-
dm <- julia_call("LarremoreModel", params$gamma, frac_symptomatic = params$frac_symptomatic, need_return = "Julia")
4+
dm <- julia_call("LarremoreModel", gamma(params$rzero, params), frac_symptomatic = params$frac_symptomatic, need_return = "Julia")
55
individuals <- julia_call("Individual.", dm, rep(params$pr_noncovid_symptoms, nrsmpl), need_return = "Julia")
66
julia_call("infect!.", individuals, need_return = "Julia")
77
for (i in 1:21) {

0 commit comments

Comments
 (0)