From 9924a3eca162f841355927932cb2d62c21e080c5 Mon Sep 17 00:00:00 2001 From: Kevin Kunzmann Date: Fri, 5 Mar 2021 16:14:01 +0000 Subject: [PATCH] refactor, 2nd preprint release --- .julia/Manifest.toml | 30 +- code/evaluate_performance.jl | 49 +- code/expand_scenarios.R | 35 -- code/scenario.R | 182 ++++++++ code/school.R | 16 + code/school.jl | 19 +- code/setup.jl | 5 +- code/simulate_r_zero.jl | 12 - code/util.R | 7 +- plots/autocorrelation.R | 65 +++ plots/calibration-infectivity.R | 36 ++ plots/fit-mean-sensitivity.R | 36 ++ plots/lower-lfd-compliance.R | 76 ++++ plots/lower-lli.R | 56 +++ plots/main-results.R | 95 ++++ plots/only-one-bubble.R | 58 +++ plots/population-structure.R | 38 ++ plots/sensitivity-vs-infectivity.R | 35 ++ plots/vl-trajectories.R | 31 ++ renv.lock | 21 + run.R | 696 +++-------------------------- 21 files changed, 847 insertions(+), 751 deletions(-) delete mode 100644 code/expand_scenarios.R create mode 100644 code/scenario.R create mode 100644 code/school.R delete mode 100644 code/simulate_r_zero.jl create mode 100644 plots/autocorrelation.R create mode 100644 plots/calibration-infectivity.R create mode 100644 plots/fit-mean-sensitivity.R create mode 100644 plots/lower-lfd-compliance.R create mode 100644 plots/lower-lli.R create mode 100644 plots/main-results.R create mode 100644 plots/only-one-bubble.R create mode 100644 plots/population-structure.R create mode 100644 plots/sensitivity-vs-infectivity.R create mode 100644 plots/vl-trajectories.R diff --git a/.julia/Manifest.toml b/.julia/Manifest.toml index b74dca6..b92dff5 100644 --- a/.julia/Manifest.toml +++ b/.julia/Manifest.toml @@ -92,9 +92,9 @@ uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b" [[Distributions]] deps = ["FillArrays", "LinearAlgebra", "PDMats", "Printf", "QuadGK", "Random", "SparseArrays", "SpecialFunctions", "Statistics", "StatsBase", "StatsFuns"] -git-tree-sha1 = "0fc424e725eaec6ea3e9fa8df773bee18a1ab503" +git-tree-sha1 = "e64debe8cd174cc52d7dd617ebc5492c6f8b698c" uuid = "31c24e10-a181-5473-b8eb-7969acd0382f" -version = "0.24.14" +version = "0.24.15" [[DocStringExtensions]] deps = ["LibGit2", "Markdown", "Pkg", "Test"] @@ -104,9 +104,9 @@ version = "0.8.3" [[Documenter]] deps = ["Base64", "Dates", "DocStringExtensions", "IOCapture", "InteractiveUtils", "JSON", "LibGit2", "Logging", "Markdown", "REPL", "Test", "Unicode"] -git-tree-sha1 = "21fb992ef1b28ff8f315354d3808ebf4a8fa6e45" +git-tree-sha1 = "3ebb967819b284dc1e3c0422229b58a40a255649" uuid = "e30172f5-a6a5-5a46-863b-614d45cd2de4" -version = "0.26.2" +version = "0.26.3" [[FillArrays]] deps = ["LinearAlgebra", "Random", "SparseArrays"] @@ -191,9 +191,9 @@ uuid = "a63ad114-7e13-5084-954f-fe012c677804" [[OffsetArrays]] deps = ["Adapt"] -git-tree-sha1 = "f64fbf703e7f5af5d82ea7b7385b443bb7765ce9" +git-tree-sha1 = "b3dfef5f2be7d7eb0e782ba9146a5271ee426e90" uuid = "6fe1bfb0-de20-5000-8ca7-80f57d26f881" -version = "1.6.1" +version = "1.6.2" [[OpenSpecFun_jll]] deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"] @@ -207,26 +207,26 @@ uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d" version = "1.4.0" [[PDMats]] -deps = ["LinearAlgebra", "SparseArrays", "SuiteSparse", "Test"] -git-tree-sha1 = "95a4038d1011dfdbde7cecd2ad0ac411e53ab1bc" +deps = ["LinearAlgebra", "SparseArrays", "SuiteSparse"] +git-tree-sha1 = "f82a0e71f222199de8e9eb9a09977bd0767d52a0" uuid = "90014a1f-27ba-587c-ab20-58faa44d9150" -version = "0.10.1" +version = "0.11.0" [[Parsers]] deps = ["Dates"] -git-tree-sha1 = "50c9a9ed8c714945e01cd53a21007ed3865ed714" +git-tree-sha1 = "223a825cccef2228f3fdbf2ecc7ca93363059073" uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0" -version = "1.0.15" +version = "1.0.16" [[Pkg]] deps = ["Dates", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "UUIDs"] uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" [[PooledArrays]] -deps = ["DataAPI"] -git-tree-sha1 = "0e8f5c428a41a81cd71f76d76f2fc3415fe5a676" +deps = ["DataAPI", "Future"] +git-tree-sha1 = "cde4ce9d6f33219465b55162811d8de8139c0414" uuid = "2dfb63ee-cc39-5dd5-95bd-886bf059d720" -version = "1.1.0" +version = "1.2.1" [[PrettyTables]] deps = ["Crayons", "Formatting", "Markdown", "Reexport", "Tables"] @@ -367,7 +367,7 @@ version = "0.5.3" [[cov19sim]] deps = ["DataFrames", "DataFramesMeta", "Distributed", "Distributions", "Documenter", "Interpolations", "Random", "Statistics", "Test", "UUIDs"] -git-tree-sha1 = "788fc38f4969992275a76e179e8bf713c412c9ed" +git-tree-sha1 = "d200e3affaab6a4b371ef354dcb94797992c56c0" repo-rev = "main" repo-url = "https://github.com/kkmann/cov19sim" uuid = "7070818b-bc4a-4efc-9908-1d2131e194b1" diff --git a/code/evaluate_performance.jl b/code/evaluate_performance.jl index 136fad0..3f816bc 100644 --- a/code/evaluate_performance.jl +++ b/code/evaluate_performance.jl @@ -1,39 +1,14 @@ -function evaluate_performance_( - n_per_bubble, bubbles_per_class, m_classes, pr_meet_class, pr_meet_school, - gamma, frac_symptomatic, pr_external_infections, pr_noncovid_symptoms, - pcr, - eta, slope, intercept, ar_window, ar_coefficient, - policy_name, a, b -) - lfd = LogRegTest("lfd", eta*slope, intercept, .998; ar_window = ar_window, ar_coefficient = ar_coefficient) - s = school( - n_per_bubble, bubbles_per_class, m_classes, pr_meet_class, pr_meet_school, - gamma, frac_symptomatic, pr_noncovid_symptoms, a, b; - policy = get_policy(policy_name, lfd) - ) - for i in 1:(6*7) - for indv in s.individuals - if rand(Bernoulli(pr_external_infections)) - infect!(indv) +function f(school, pr_external_infections, days; n = 1L) + function g(x) + for i in 1:(days) + for indv in x.individuals + if rand(Bernoulli(pr_external_infections)) + infect!(indv) + end end - end - step!(s) - end - evaluate(s) + step!(x) + end + return evaluate(x) + end + return vcat(pmap(g, [resample(school) for i = 1:n])...) end - -function evaluate_performance( - n_per_bubble, bubbles_per_class, m_classes, pr_meet_class, pr_meet_school, - gamma, frac_symptomatic, pr_external_infections, pr_noncovid_symptoms, - pcr, - eta, slope, intercept, ar_window, ar_coefficient, - policy_name, a, b -) - vcat(pmap(evaluate_performance_, - n_per_bubble, bubbles_per_class, m_classes, pr_meet_class, pr_meet_school, - gamma, frac_symptomatic, pr_external_infections, pr_noncovid_symptoms, - pcr, - eta, slope, intercept, ar_window, ar_coefficient, - policy_name, a, b - )...) -end \ No newline at end of file diff --git a/code/expand_scenarios.R b/code/expand_scenarios.R deleted file mode 100644 index 4dfdbdc..0000000 --- a/code/expand_scenarios.R +++ /dev/null @@ -1,35 +0,0 @@ -expand_scenarios <- function( - n_bubble = 9, - bubbles_per_class = 3, - classes_per_school = 12, - pr_meet_class = 3/(n_bubble*bubbles_per_class - 1), - pr_meet_school = 1/(n_bubble*bubbles_per_class*classes_per_school - 1), - R = c(1.5, 3, 6), - frac_symptomatic = 0.5, - pr_external_infections = 1/(n_bubble*bubbles_per_class*classes_per_school*7), # expected 1 per week - pr_noncovid_symptoms = .01, - mean_sensitivity = c(.4, .6, .8), - ar_coefficient = 0.0, - ar_window = 3L, - a = Inf, - b = 1.0, - pcr_test = list(pcr), - slope = innova$slope, - intercept = innova$intercept, - iter = 1:(params$iterations), - policy_name = c("default", "Thu/Fri off", "Mon screening", "Mon/Wed screening", "test & release") -) { - args <- as.list(environment()) - args[["R"]] <- tibble(R = R, gamma = purrr::map_dbl(R, get_gamma)) - args[["mean_sensitivity"]] <- tibble( - mean_sensitivity = mean_sensitivity, - eta = purrr::map_dbl(mean_sensitivity, get_eta) - ) - do.call(expand_grid, args) %>% - mutate( - gamma = R$gamma, - R = R$R, - eta = mean_sensitivity$eta, - mean_sensitivity = mean_sensitivity$mean_sensitivity - ) -} diff --git a/code/scenario.R b/code/scenario.R new file mode 100644 index 0000000..6ab4ab8 --- /dev/null +++ b/code/scenario.R @@ -0,0 +1,182 @@ +scenario <- function( + n_bubble = 9L, + bubbles_per_class = 3L, + classes_per_school = 12L, + expected_class_contacts = 3, + expected_school_contacts = 1, + frac_symptomatic = 0.5, + pr_noncovid_symptoms = 0.01, + expected_weekly_external_infections = 1, + a = Inf, + b = 1, + lli = 1e6, + pcr_lod = 300.0, + pcr_sens = .975, + pcr_spec = 1.0, + gamma = 0.0, + eta = 1.0, + lfd_spec = 0.998, + ar_window = 3L, + ar_coefficient = 0.0, + days = as.integer(6*7), + ... +) { + res <- as.list(environment()) + + res$pr_meet_class <- expected_class_contacts/(n_class(res) - 1) + res$expected_class_contacts <- NULL + res$pr_meet_school <- expected_school_contacts/(n_school(res) - 1) + res$expected_school_contacts <- NULL + res$pr_external_infection <- max(0, min(1, expected_weekly_external_infections / n_school(res) / 7)) + res$expected_weekly_external_infections <- NULL + + tbl_innova_data <- tibble( + `viral load` = 10^(2:7), + sensitivity = 1 - (c(621, 521, 344, 207, 144, 123) - 114) / (684 - 114) + ) + sensitivity <- function(vl, slope, intercept) { + lfd <- julia_call("LogRegTest", "lfd", slope, intercept, lfd_spec, need_return = "Julia") + julia_call("sensitivity.", lfd, vl, need_return = "R") + } + innova <- optim( + c(.76, -5), + function(x) { + sum(( + tbl_innova_data$sensitivity - + sensitivity(tbl_innova_data$`viral load`, slope = x[1], intercept = x[2]) + )^2) + }, + method = "L-BFGS-B", + lower = c(.1, -10), + upper = c(2.5, 5), + control = list(maxit = 1e4) + ) + res$lfd_slope <- innova$par[1] + res$lfd_intercept <- innova$par[2] + + res +} + +sensitivity <- function(vl, params) { + lfd <- julia_call("LogRegTest", "lfd", params$lfd_slope, params$lfd_intercept, params$lfd_spec, need_return = "Julia") + julia_call("sensitivity.", lfd, vl, need_return = "R") +} + +n_class <- function(params) with(params, n_bubble * bubbles_per_class) +n_school <- function(params) with(params, classes_per_school * n_class(params)) +n_weekly_infections <- function(params) with(params, 7*n_school(params)*pr_external_infection) + +schooldays <- function(params) with(params, sum(((1:days) %% 7) %in% 1:5)) + +pcr <- function(params) with(params, julia_call("FixedTest", "pcr", pcr_sens, pcr_spec, lod = pcr_lod, need_return = "Julia")) + +lfd <- function(params) { + with(params, + julia_call("LogRegTest", "lfd", + eta*params$lfd_slope, params$lfd_intercept, lfd_spec, + ar_window = ar_window, ar_coefficient = ar_coefficient + ) + ) +} + +fit_rzero_ <- function(params, gamma_min = 0, gamma_max = .1) { + tbl_rzero <- tibble( + gamma = seq(gamma_min, gamma_max, length.out = 100), + R = map(gamma, ~{ + params$gamma <- . + sample_rzero(do.call(school, params), n = 10L) + }) + ) %>% + unnest(R) + fit <- lm(formula = R ~ gamma - 1, data = tbl_rzero) + list(fit = fit, data = tbl_rzero) +} +fit_rzero <- memoise::memoise(fit_rzero_) +rzero <- function(gamma, params, gamma_min = 0, gamma_max = .1) { + fit <- fit_rzero(params, gamma_min, gamma_max)$fit + as.numeric(predict(fit, newdata = tibble(gamma = gamma), type = "response")) +} +gamma <- function(R, params, gamma_min = 0, gamma_max = .1) { + as.numeric(uniroot(function(x) rzero(x, params, gamma_min, gamma_max) - R, interval = c(gamma_min, gamma_max))$root) +} + +sample_presymptomatic_vl <- function(params, n = 1e5) { + dm <- julia_call("LarremoreModel", params$gamma, frac_symptomatic = params$frac_symptomatic, need_return = "Julia") + individuals <- julia_call("Individual.", dm, rep(params$pr_noncovid_symptoms, n), need_return = "Julia") + julia_call("infect!.", individuals, need_return = "Julia") + julia_call("steps!.", individuals, 21L, need_return = "Julia") + julia_call("get_status_logs", individuals) %>% + as_tibble() %>% + arrange(uuid, day) %>% + group_by(uuid) %>% + filter( + row_number() >= which(viral_load > params$pcr_lod)[1], + row_number() < which(symptomatic)[1] + ) %>% + sample_n(1) %>% + ungroup() %>% + select(uuid, day, viral_load) +} +sample_presymptomatic_vl_ <- memoise::memoise(sample_presymptomatic_vl) +mean_sensitivity <- function(params, eta, ...) { + sample_presymptomatic_vl_(params, ...) %>% + mutate( + sensitivity = sensitivity(viral_load^eta, params) + ) %>% + pull(sensitivity) %>% + summary() %>% + mean() +} +eta <- function(params, target, ...) { + as.numeric(uniroot(function(x) mean_sensitivity(params, x, ...) - target, interval = c(0, 5))$root) +} + +expand_scenario <- function(params = scenario(), ...) { + tbl <- if (length(list(...)) == 0) { + as_tibble(params) + } else { + for (name in names(list(...))) { + params[[name]] <- NULL + } + expand_grid(do.call(expand_grid, list(...)), as_tibble(params)) + } + tbl %>% + group_by(row_number()) %>% + nest() %>% + mutate(data = map(data, as.list)) %>% + pull(data) %>% + enframe(value = "params") %>% + select(params) %>% + mutate( + tmp = map(params, ~as_tibble(.)) + ) %>% + unnest(tmp) +} + +julia_eval('@everywhere include("code/evaluate_performance.jl")') +evaluate_performance_ <- function(params, policy, n = 1L) { + school <- do.call(school, c(params, list(policy = policy(params)))) + julia_call("f", school, params$pr_external_infection, params$days, n = n) %>% + as_tibble() %>% + mutate( + `% infected (cumulative)` = n_infected/n_school(params), + `% schooldays missed (cumulative)` = workdays_missed/n_school(params)/schooldays(params) + ) +} +evaluate_performance <- function(policies, params = scenario(), n = if (!is.null(n_resample)) {n_resample} else {25L}, ...) { + gamma2rs <- memoise::memoise(function(gamma) round(rzero(gamma, params), 1) ) + eta2mean_sensitivity <- memoise::memoise(function(eta) round(mean_sensitivity(params, eta), 2) ) + expand_scenario(params, ...) %>% + expand_grid(tibble(policy = policies)) %>% + mutate( + policy_name = names(policy), + results = map2(params, policy, evaluate_performance_, n) + ) %>% + unnest(results) %>% + select(-params) %>% + mutate( + policy_name = factor(policy_name, levels = names(lst_policies)), + Rs = map_dbl(.data$gamma, gamma2rs), + `mean sensitivity` = map_dbl(.data$eta, eta2mean_sensitivity) + ) +} diff --git a/code/school.R b/code/school.R new file mode 100644 index 0000000..e5fa8f3 --- /dev/null +++ b/code/school.R @@ -0,0 +1,16 @@ +julia_eval('@everywhere include("code/school.jl")') + +school <- function(policy = julia_call("DoNothing", need_return = "Julia"), ...){ + params <- scenario(...) + with(params, + julia_call("school", + n_bubble, bubbles_per_class, classes_per_school, pr_meet_class, pr_meet_school, + gamma, + frac_symptomatic, pr_noncovid_symptoms, + a, b, + lli, + policy + )) +} + +sample_rzero <- function(school, days = 21L, n = 1L) julia_call("sample_rzero", school, days, n) diff --git a/code/school.jl b/code/school.jl index fd86297..08ed4e1 100644 --- a/code/school.jl +++ b/code/school.jl @@ -1,9 +1,9 @@ function school( n_per_bubble, bubbles_per_class, m_classes, pr_meet_class, pr_meet_school, - gamma, frac_symptomatic, pr_noncov_symptoms, a, b; - policy = DoNothing() + gamma, frac_symptomatic, pr_noncov_symptoms, a, b, lli, + policy ) - dm = LarremoreModel(gamma; frac_symptomatic = frac_symptomatic) + dm = LarremoreModel(gamma; frac_symptomatic = frac_symptomatic, l10vl_clearance = log10(lli)) ThreeLevelPopulation( n_per_bubble = n_per_bubble, bubbles_per_class = bubbles_per_class, @@ -20,4 +20,17 @@ function school( a = a, b = b ) +end + +function sample_rzero(school::T, days::Int, n::Int) where {T<:Population} + function f() + m = length(school.individuals) + tmp = resample(school) + ind = rand(1:m) + infect!.(tmp.individuals[ind]) + steps!(tmp, days) + res = get_contact_log(tmp.individuals[ind]) + return nrow(res[res[:, :infected_other], :]) + end + pmap(i -> f(), 1:n) end \ No newline at end of file diff --git a/code/setup.jl b/code/setup.jl index cb13540..6498e91 100644 --- a/code/setup.jl +++ b/code/setup.jl @@ -10,14 +10,13 @@ seed!(42) # load environment on worker processors using Distributed -addprocs(min(8, Sys.CPU_THREADS)) # memory limitations... +rmprocs(workers()) +addprocs(Sys.CPU_THREADS) # memory limitations... @everywhere begin using Pkg Pkg.activate(".julia") using cov19sim, DataFrames, Distributions, Random Random.seed!(myid()) include("code/school.jl") - include("code/simulate_r_zero.jl") - include("code/get_policy.jl") include("code/evaluate_performance.jl") end \ No newline at end of file diff --git a/code/simulate_r_zero.jl b/code/simulate_r_zero.jl deleted file mode 100644 index 28d7f82..0000000 --- a/code/simulate_r_zero.jl +++ /dev/null @@ -1,12 +0,0 @@ -function simulate_r_zero(school::T; days = 21) where {T<:Population} - school = resample(school) - infect!.(school.individuals[1]) - steps!(school, days) - res = get_contact_log(school.individuals[1]) - nrow(res[res[:, :infected_other], :]) -end - -# pmap vectorized form -function simulate_r_zero(schools::Vector{T}; days = 21) where {T<:Population} - pmap(x -> simulate_r_zero(x, days = days), schools) -end \ No newline at end of file diff --git a/code/util.R b/code/util.R index 832c043..83a78a4 100644 --- a/code/util.R +++ b/code/util.R @@ -2,10 +2,6 @@ save_plot <- function(plt, name, width = 7.5, height = width*3/4, dpi = 300, ... dir <- "_site/figures" if (!dir.exists(dir)) dir.create(dir, recursive = TRUE) - plt <- plt + ggplot2::labs(...) # allow removal of caption etc. - ggplot2::ggsave(sprintf("%s/%s.jpg", dir, name), plt, - width = width, height = height, dpi = dpi) - plt <- plt + ggplot2::labs(caption = NULL) # no caption for pdf ggplot2::ggsave(sprintf("%s/%s.pdf", dir, name), plt, width = width, height = height, dpi = dpi) } @@ -16,6 +12,7 @@ theme_set( legend.position = "top", panel.grid.minor = element_blank(), legend.title = element_blank(), - plot.caption = element_text(size = 11, hjust = .5, vjust = 0) + plot.caption = element_text(size = 11, hjust = .5, vjust = 0), + panel.spacing = unit(1.25, "lines") ) ) diff --git a/plots/autocorrelation.R b/plots/autocorrelation.R new file mode 100644 index 0000000..b87a64a --- /dev/null +++ b/plots/autocorrelation.R @@ -0,0 +1,65 @@ +get_ar <- function(params, seed = 42L, nrsmpl = 1e4) { + julia_call("Random.seed!", as.integer(seed)) + lfd <- lfd(params) + dm <- julia_call("LarremoreModel", params$gamma, frac_symptomatic = params$frac_symptomatic, need_return = "Julia") + individuals <- julia_call("Individual.", dm, rep(params$pr_noncovid_symptoms, nrsmpl), need_return = "Julia") + julia_call("infect!.", individuals, need_return = "Julia") + for (i in 1:21) { + if ((i %% 7) %in% 0:4) { # only test on weekends + julia_call("conduct_test!.", lfd, individuals, need_return = "None") + } + julia_call("step!.", individuals, need_return = "None") + } + tbl_status <- julia_call("get_status_logs", individuals, need_return = "R") %>% + as_tibble() + tbl_tests <- julia_call("get_test_logs", individuals, need_return = "R") %>% + as_tibble() %>% + select(-viral_load) + + left_join(tbl_tests, tbl_status, by = c("uuid", "day")) %>% + select(uuid, day, result, viral_load, probability_positive, symptomatic) %>% + group_by(uuid) %>% + filter( + !any(symptomatic) | (row_number() < which(symptomatic)[1]) + ) %>% + summarize( + acf = list(tibble( + lag = acf(as.numeric(result), plot = FALSE, lag.max = 7)$lag[,,1], + ac = acf(as.numeric(result), plot = FALSE, lag.max = 7)$acf[,,1] + )), + .groups = "drop" + ) %>% + unnest(acf) %>% + mutate( + ac = if_else(is.nan(ac), 1, ac) # constant signal has ac 1 + ) %>% + group_by(lag) %>% + summarize( + mid = mean(ac), + lo = quantile(ac, .05), + hi = quantile(ac, .95), + .groups = "drop" + ) +} + +plt <- tibble( + coefficient = c(0, .25, .5, .75) + ) %>% + mutate( + ar = map(coefficient, ~get_ar(scenario(ar_coefficient = .), nrsmpl = 1e4)) + ) %>% + unnest(ar) %>% + mutate( + coefficient = sprintf("%.2f", coefficient) + ) %>% + ggplot() + + aes(lag, mid, color = coefficient) + + geom_line(aes(y = mid)) + + scale_color_grey("") + + ylab("autocorrelation") + + scale_x_continuous("lag [days]", breaks = seq(0, 21, by = 1)) + + theme( + legend.position = "right" + ) + +save_plot(plt, "retest-autocorrelation", width = width, height = height/2) \ No newline at end of file diff --git a/plots/calibration-infectivity.R b/plots/calibration-infectivity.R new file mode 100644 index 0000000..8338443 --- /dev/null +++ b/plots/calibration-infectivity.R @@ -0,0 +1,36 @@ +params1 <- scenario() +plt1 <- ggplot(fit_rzero(params1)$data) + + aes(gamma, R) + + geom_point(alpha = 0.2, shape = 16) + + geom_line( + data = tibble( + gamma = seq(0, .1, length.out = 100), + R = as.numeric(predict( fit_rzero(params1)$fit, newdata = tibble(gamma = gamma), type = "response")) + ) + ) + + ylim(c(0, 16)) + + xlim(c(0, .1)) + + labs( + y = "", + x = expression(gamma), + title = "LLI = 1e6" + ) +params2 <- scenario(lli = 1e3) +plt2 <- ggplot(fit_rzero(params2, gamma_max = .06)$data) + + aes(gamma, R) + + geom_point(alpha = 0.2, shape = 16) + + geom_line( + data = tibble( + gamma = seq(0, .06, length.out = 100), + R = as.numeric(predict( fit_rzero(params2, gamma_max = .06)$fit, newdata = tibble(gamma = gamma), type = "response")) + ) + ) + + ylim(c(0, 16)) + + xlim(c(0, .1)) + + labs( + y = "", + x = expression(gamma), + title = "LLI = 1e3" + ) +plt <- plt1 + plt2 +save_plot(plt, "calibrating-infectivity", width = width, height = .66*height) diff --git a/plots/fit-mean-sensitivity.R b/plots/fit-mean-sensitivity.R new file mode 100644 index 0000000..8b6e079 --- /dev/null +++ b/plots/fit-mean-sensitivity.R @@ -0,0 +1,36 @@ +params <- scenario() + +tbl_innova_data <- tibble( + `viral load` = 10^(2:7), + sensitivity = 1 - (c(621, 521, 344, 207, 144, 123) - 114) / (684 - 114) + ) +tbl_scenarios <- tibble( + `mean sensitivity` = c(mean_sensitivity(scenario(), 1), .4, .6, .8), + eta = map_dbl(`mean sensitivity`, ~eta(scenario(), target = .)) + ) +plt <- tbl_scenarios %>% + expand_grid( + `viral load` = 10^seq(0, 11, length.out = 1000), + ) %>% + mutate( + sensitivity = sensitivity(`viral load`^eta, params) + ) %>% + ggplot() + + aes(`viral load`, sensitivity) + + geom_line(aes(linetype = sprintf("%.2f", `mean sensitivity`))) + + geom_point(data = tbl_innova_data) + + scale_x_log10() + + scale_y_continuous("", limits = c(0, 1)) + + scale_linetype_manual("", + breaks = sprintf("%.2f", tbl_scenarios$`mean sensitivity`), + values = c(1, 4, 3, 2), + labels = purrr::map2( + sprintf("%.2f", tbl_scenarios$`mean sensitivity`), + tbl_scenarios$eta, + ~bquote( .(..1) (eta == .(round(..2, 2))) ) + ) + ) + + theme( + legend.position = "right" + ) +save_plot(plt, "fit-mean-sensitivity", width = width, height = height/2) \ No newline at end of file diff --git a/plots/lower-lfd-compliance.R b/plots/lower-lfd-compliance.R new file mode 100644 index 0000000..e9827c8 --- /dev/null +++ b/plots/lower-lfd-compliance.R @@ -0,0 +1,76 @@ +plt <- tibble( + compliance = rbeta(10^5, shape1 = 2/15, shape2 = 1/15) + ) %>% + ggplot() + + aes(compliance) + + geom_histogram(aes(y = ..ncount..), binwidth = 0.05) + + scale_y_continuous("", labels = scales::percent) +save_plot(plt, "beta-distribution", width = width, height = .66*height) + + +params <- scenario(a = 2/15, b = 1/15) +tbl_results <- bind_rows( + evaluate_performance( + policies = lst_policies[c("Mon/Wed screening", "Mon screening", "test for release")], + params = scenario(), # standard compliance + gamma = map_dbl(c(1.5, 3, 6), ~gamma(., scenario())), + eta = map_dbl(c(0.4, 0.6, 0.8), ~eta(scenario(), .)) + ), + evaluate_performance( + policies = lst_policies[c("Mon/Wed screening", "Mon screening", "test for release")], + params = params, + gamma = map_dbl(c(1.5, 3, 6), ~gamma(., params)), + eta = map_dbl(c(0.4, 0.6, 0.8), ~eta(params, .)) + ) +) +dir.create("_site/data", showWarnings = FALSE, recursive = TRUE) +write_rds(tbl_results, "_site/data/tbl_compliance.rds") + +plt1 <- tbl_results %>% + mutate( + tmp = if_else(is.finite(a), a / (a + b), 1), + `mean LFD test compliance:` = factor( + tmp, + levels = unique(tmp), + labels = sprintf("%.2f", unique(tmp)) + ) + ) %>% + ggplot() + + aes(policy_name, `% infected (cumulative)`, color = `mean LFD test compliance:`) + + geom_boxplot() + + scale_y_continuous(labels = scales::percent, limits = c(0, NA_real_)) + + facet_wrap(~Rs, labeller = label_both, nrow = 1) + + theme( + axis.text.x = element_text(angle = 33, hjust = 1), + axis.title.x = element_blank(), + legend.position = "top", + legend.title = element_text() + ) + +plt2 <- tbl_results %>% + mutate( + tmp = if_else(is.finite(a), a / (a + b), 1), + `mean LFD test compliance:` = factor( + tmp, + levels = unique(tmp), + labels = sprintf("%.2f", unique(tmp)) + ) + ) %>% + filter( + `mean LFD test compliance:` == "0.67" + ) %>% + ggplot() + + aes(`% infected (cumulative)`, `% schooldays missed (cumulative)`, color = policy_name, fill = policy_name) + + facet_grid(`mean sensitivity` ~ Rs, labeller = label_both) + + geom_abline(slope = 1) + + geom_point(alpha = .33, shape = 16, size = 2) + + scale_x_continuous(labels = scales::percent, limits = c(0, NA_real_)) + + scale_y_continuous(labels = scales::percent, limits = c(0, NA_real_)) + + scale_color_discrete(guide = guide_legend(override.aes = list(alpha = 1))) + + coord_cartesian(expand = FALSE) + + theme( + legend.text = element_text(size = rel(1.1)) + ) +plt <- plt1 + plt2 + plot_layout(ncol = 1, heights = c(1, 3)) + plot_annotation(tag_levels = "A") +save_plot(plt, "results-schooldays-missed-vs-infectivity-lower-lfd-compliance", width = width, height = 1.75*height) + diff --git a/plots/lower-lli.R b/plots/lower-lli.R new file mode 100644 index 0000000..f585bc0 --- /dev/null +++ b/plots/lower-lli.R @@ -0,0 +1,56 @@ +params <- scenario(lli = 1e3) +tbl_results <- bind_rows( + evaluate_performance( + policies = lst_policies, + params = scenario(), + gamma = map_dbl(c(1.5, 3, 6), ~gamma(., scenario())), + eta = eta(scenario(), target = 0.6) + ), + evaluate_performance( + policies = lst_policies, + params = params, + gamma = map_dbl(c(1.5, 3, 6), ~gamma(., params, gamma_max = .06)), + eta = eta(params, target = 0.6) + ) %>% + mutate( + Rs = map_dbl(.data$gamma, ~round(rzero(., params, gamma_max = .06), 1)) + ) + ) %>% + mutate( + LLI = factor(lli, levels = c(scenario()$lli, params$lli)) + ) +dir.create("_site/data", showWarnings = FALSE, recursive = TRUE) +write_rds(tbl_results, "_site/data/tbl_lli.rds") + +plt1 <- tbl_results %>% + ggplot() + + aes(policy_name, `% infected (cumulative)`, color = LLI) + + geom_boxplot() + + scale_y_continuous(labels = scales::percent, limits = c(0, NA_real_)) + + scale_color_discrete("LLI:") + + facet_wrap(~Rs, labeller = label_both, nrow = 1) + + theme( + axis.text.x = element_text(angle = 33, hjust = 1), + axis.title.x = element_blank(), + legend.position = "top", + legend.title = element_text() + ) + +plt2 <- tbl_results %>% + filter( + LLI == params$lli + ) %>% + ggplot() + + aes(`% infected (cumulative)`, `% schooldays missed (cumulative)`, color = policy_name, fill = policy_name) + + facet_grid(LLI ~ Rs, labeller = label_both) + + geom_abline(slope = 1) + + geom_point(alpha = .33, shape = 16, size = 2) + + scale_x_continuous(labels = scales::percent, limits = c(0, NA_real_)) + + scale_y_continuous(labels = scales::percent, limits = c(0, NA_real_)) + + scale_color_discrete(guide = guide_legend(override.aes = list(alpha = 1))) + + coord_cartesian(expand = FALSE) + + theme( + legend.text = element_text(size = rel(1.1)) + ) +plt <- plt1 + plt2 + plot_layout(ncol = 1, heights = c(1, 1.5)) + plot_annotation(tag_levels = "A") +save_plot(plt, "results-schooldays-missed-vs-infectivity-lower-lli", width = width, height = 1.5*height) diff --git a/plots/main-results.R b/plots/main-results.R new file mode 100644 index 0000000..36bf3bd --- /dev/null +++ b/plots/main-results.R @@ -0,0 +1,95 @@ +params <- scenario() +tbl_results <- evaluate_performance( + policies = lst_policies, + params = params, + gamma = map_dbl(c(1.5, 3, 6), ~gamma(., params)), + eta = eta(params, target = 0.6), + frac_symptomatic = c(0.25, 0.5, 0.75) + ) +dir.create("_site/data", showWarnings = FALSE, recursive = TRUE) +write_rds(tbl_results, "_site/data/tbl_main.rds") + +plt <- tbl_results %>% + filter( + Rs == 3.0, + frac_symptomatic == 0.5 + ) %>% + select( + policy_name, + `% infected (cumulative)`, + `% schooldays missed (cumulative)`, + ) %>% + pivot_longer(-policy_name) %>% + ggplot() + + aes(policy_name, value) + + geom_boxplot() + + scale_y_continuous("", labels = scales::percent, limits = c(0, NA_real_)) + + facet_wrap(~name) + + theme( + axis.text.x = element_text(angle = 33, hjust = 1), + axis.title.x = element_blank(), + legend.position = "right" + ) +save_plot(plt, "results-main", width = width, height = .75*height) + +plt1 <- tbl_results %>% + ggplot() + + aes(policy_name, `% infected (cumulative)`, color = factor(1 - frac_symptomatic)) + + geom_boxplot() + + scale_y_continuous(labels = scales::percent, limits = c(0, NA_real_)) + + scale_color_discrete("% asymptomatic cases:") + + facet_grid(`mean sensitivity` ~ Rs, labeller = label_both) + + theme( + axis.text.x = element_text(angle = 33, hjust = 1), + axis.title.x = element_blank(), + legend.position = "top", + legend.title = element_text() + ) + + +tbl_results2 <- evaluate_performance( + policies = lst_policies[c("Mon/Wed screening", "Mon screening", "test for release")], + params = params, + gamma = map_dbl(c(3), ~gamma(., params)), + eta = map_dbl(c(.4, .6, .8), ~eta(params, .)), + ar_coefficient = c(0, 0.75) + ) +dir.create("_site/data", showWarnings = FALSE, recursive = TRUE) +write_rds(tbl_results, "_site/data/tbl_ar_coefficient.rds") +plt2 <- tbl_results2 %>% + ggplot() + + aes(policy_name, `% infected (cumulative)`, color = factor(ar_coefficient)) + + geom_boxplot() + + scale_y_continuous(labels = scales::percent, limits = c(0, NA_real_)) + + scale_color_discrete("AR coefficient:") + + facet_grid(Rs ~ `mean sensitivity`, labeller = label_both) + + theme( + axis.text.x = element_text(angle = 33, hjust = 1), + axis.title.x = element_blank(), + legend.position = "top", + legend.title = element_text() + ) + +plt <- plt1 + plt2 + plot_layout(nrow = 2) + plot_annotation(tag_levels = "A") +save_plot(plt, "results-symptoms-ar", width = width, height = 1.25*height) + + + + +plt <- tbl_results %>% + filter( + frac_symptomatic == .5 + ) %>% + ggplot() + + aes(`% infected (cumulative)`, `% schooldays missed (cumulative)`, color = policy_name, fill = policy_name) + + facet_grid(`mean sensitivity` ~ Rs, labeller = label_both) + + geom_abline(slope = 1) + + geom_point(alpha = .33, shape = 16, size = 2) + + scale_x_continuous(labels = scales::percent, limits = c(0, NA_real_)) + + scale_y_continuous(labels = scales::percent, limits = c(0, NA_real_)) + + scale_color_discrete(guide = guide_legend(override.aes = list(alpha = 1))) + + coord_cartesian(expand = FALSE) + + theme( + legend.text = element_text(size = rel(1.1)) + ) +save_plot(plt, "results-schooldays-missed-vs-infectivity", width = width, height = .75*height) diff --git a/plots/only-one-bubble.R b/plots/only-one-bubble.R new file mode 100644 index 0000000..a33246c --- /dev/null +++ b/plots/only-one-bubble.R @@ -0,0 +1,58 @@ +params <- scenario(n_bubble = 27, bubbles_per_class = 1) +tbl_results <- bind_rows( + evaluate_performance( + policies = lst_policies, + params = scenario(), + gamma = map_dbl(c(1.5, 3, 6), ~gamma(., scenario())), + eta = eta(scenario(), target = 0.6) + ), + evaluate_performance( + policies = lst_policies, + params = params, + gamma = map_dbl(c(1.5, 3, 6), ~gamma(., scenario())), + eta = eta(scenario(), target = 0.6) + ) +) +dir.create("_site/data", showWarnings = FALSE, recursive = TRUE) +write_rds(tbl_results, "_site/data/tbl_one_bubble.rds") + +gamma2rs <- memoise::memoise(function(gamma) round(rzero(gamma, scenario()), 1)) + +plt1 <- tbl_results %>% + mutate( + `bubbles per class` = factor(bubbles_per_class, levels = c(3, 1)), + Rs = map_dbl(gamma, gamma2rs) + ) %>% + ggplot() + + aes(policy_name, `% infected (cumulative)`, color = `bubbles per class`) + + geom_boxplot() + + scale_y_continuous(labels = scales::percent, limits = c(0, NA_real_)) + + facet_wrap(~Rs, labeller = label_both, nrow = 1) + + theme( + axis.text.x = element_text(angle = 33, hjust = 1), + axis.title.x = element_blank(), + legend.position = "top", + legend.title = element_text() + ) +plt2 <- tbl_results %>% + mutate( + `bubbles per class` = factor(bubbles_per_class, levels = c(3, 1)), + Rs = map_dbl(gamma, gamma2rs) + ) %>% + filter( + `bubbles per class` == 1 + ) %>% + ggplot() + + aes(`% infected (cumulative)`, `% schooldays missed (cumulative)`, color = policy_name, fill = policy_name) + + facet_grid(`bubbles per class` ~ Rs, labeller = label_both) + + geom_abline(slope = 1) + + geom_point(alpha = .33, shape = 16, size = 2) + + scale_x_continuous(labels = scales::percent, limits = c(0, NA_real_)) + + scale_y_continuous(labels = scales::percent, limits = c(0, NA_real_)) + + scale_color_discrete(guide = guide_legend(override.aes = list(alpha = 1))) + + coord_cartesian(expand = FALSE) + + theme( + legend.text = element_text(size = rel(1.1)) + ) +plt <- plt1 + plt2 + plot_layout(ncol = 1, heights = c(1, 1.5)) + plot_annotation(tag_levels = "A") +save_plot(plt, "results-schooldays-missed-vs-infectivity-only-one-bubble", width = width, height = 1.5*height) diff --git a/plots/population-structure.R b/plots/population-structure.R new file mode 100644 index 0000000..78e8985 --- /dev/null +++ b/plots/population-structure.R @@ -0,0 +1,38 @@ +plt1 <- do.call(school, scenario()) %>% + julia_call("get_adjacency_matrix", .) %>% + as_tbl_graph() %>% + mutate( + tmp = "a" + ) %>% + ggraph("matrix") + + geom_edge_tile(aes(fill = weight)) + + scale_edge_fill_continuous( + "", + low = "#dddddd", high = "#000000", na.value = "white", + limits = c(0, NA_real_) + ) + + coord_fixed(expand = FALSE) + + theme( + legend.position = "top" + ) + + labs(title = "3 bubbles per class") +plt2 <- do.call(school, scenario(n_bubble = 27, bubbles_per_class = 1)) %>% + julia_call("get_adjacency_matrix", .) %>% + as_tbl_graph() %>% + mutate( + tmp = "a" + ) %>% + ggraph("matrix") + + geom_edge_tile(aes(fill = weight)) + + scale_edge_fill_continuous( + "", + low = "#dddddd", high = "#000000", na.value = "white", + limits = c(0, NA_real_) + ) + + coord_fixed(expand = FALSE) + + theme( + legend.position = "top" + ) + + labs(title = "one bubble per class") +plt <- plt1 + plt2 + plot_layout(guides = "collect") & theme(legend.position = "right") +save_plot(plt, "adjacency-matrices", width = width, height = .66*height) diff --git a/plots/sensitivity-vs-infectivity.R b/plots/sensitivity-vs-infectivity.R new file mode 100644 index 0000000..e1a5a7c --- /dev/null +++ b/plots/sensitivity-vs-infectivity.R @@ -0,0 +1,35 @@ +params <- scenario() + +plt <- expand_grid( + tibble( + `mean sensitivity` = c(.4, .6, .8), + eta = map_dbl(`mean sensitivity`, ~eta(params, .)) + ), + tibble( + Rs = c(1.5, 3, 6), + gamma = map_dbl(Rs, ~gamma(., params)) + ) + ) %>% + expand_grid( + `viral load` = 10^seq(0, 11, length.out = 1000) + ) %>% + mutate( + sensitivity = sensitivity(`viral load`^eta, params), + `1e6` = pmax(0, pmin(1, gamma * (log10(`viral load`) - 6))), + `1e3` = pmax(0, pmin(1, gamma * (log10(`viral load`) - 3))) + ) %>% + pivot_longer(c(`1e3`, `1e6`), names_to = "LLI", values_to = "probability to infect") %>% + ggplot() + + aes(`probability to infect`, color = LLI) + + geom_line(aes(y = sensitivity)) + + facet_grid( + `mean sensitivity` ~ Rs, + labeller = label_both, scales = "free" + ) + + scale_color_discrete("LLI: ") + + labs( + y = "scaled LFD sensitivity", + x = "probability to infect during risk contact" + ) + + theme(legend.title = element_text()) +save_plot(plt, "sensitivity-vs-infectivity", width = width, height = 1.5*height) diff --git a/plots/vl-trajectories.R b/plots/vl-trajectories.R new file mode 100644 index 0000000..31884a3 --- /dev/null +++ b/plots/vl-trajectories.R @@ -0,0 +1,31 @@ +params <- scenario() + +dm <- julia_call("LarremoreModel", + gamma(R = 3, params), frac_symptomatic = params$frac_symptomatic, + l10vl_clearance = log10(params$lli), + need_return = "Julia" +) +individuals <- julia_call("Individual.", dm, rep(params$pr_noncovid_symptoms, 32), need_return = "Julia") +julia_call("infect!.", individuals, need_return = "Julia") +julia_call("steps!.", individuals, 21L, need_return = "Julia") + +plt <- julia_call("get_status_logs", individuals, need_return = "R") %>% + as_tibble() %>% + mutate( + symptomatic = factor(symptomatic, levels = c(FALSE, TRUE), labels = c("non-symptomatic", "symptomatic")) + ) %>% + ggplot() + + aes(day, viral_load) + + geom_hline(yintercept = params$lli) + + geom_bar(aes(fill = symptomatic), stat = "identity", alpha = .66, width = 1) + + scale_y_log10() + + scale_x_continuous("day post infection", breaks = seq(0, 35, by = 7)) + + scale_fill_manual(breaks = c("non-symptomatic", "symptomatic"), values = c("darkgray", "black")) + + ylab("viral load") + + facet_wrap(~uuid, nrow = 4) + + theme( + strip.text = element_blank(), + strip.background = element_blank(), + panel.spacing = unit(.25, "lines") + ) +save_plot(plt, "vl-trajectories", height = .8*height) \ No newline at end of file diff --git a/renv.lock b/renv.lock index c47696c..eaa2d26 100644 --- a/renv.lock +++ b/renv.lock @@ -128,6 +128,13 @@ "Repository": "CRAN", "Hash": "3cca8477943a8e840c3eb64a25fef9d2" }, + "cachem": { + "Package": "cachem", + "Version": "1.0.1", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "3259868ce0a3fab9456de8635cc34728" + }, "callr": { "Package": "callr", "Version": "3.5.1", @@ -247,6 +254,13 @@ "Repository": "CRAN", "Hash": "dad6793a5a1f73c8e91f1a1e3e834b05" }, + "fastmap": { + "Package": "fastmap", + "Version": "1.1.0", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "77bd60a6157420d4ffa93b27cf6a58b8" + }, "forcats": { "Package": "forcats", "Version": "0.5.1", @@ -429,6 +443,13 @@ "Repository": "CRAN", "Hash": "61e4a10781dd00d7d81dd06ca9b94e95" }, + "memoise": { + "Package": "memoise", + "Version": "2.0.0", + "Source": "Repository", + "Repository": "CRAN", + "Hash": "a0bc51650201a56d00a4798523cc91b3" + }, "mgcv": { "Package": "mgcv", "Version": "1.8-33", diff --git a/run.R b/run.R index 69861e1..e8a1033 100644 --- a/run.R +++ b/run.R @@ -1,8 +1,3 @@ -params <- list( - iterations = 100L -) - -## ----r-setup, include=FALSE, cache=FALSE---------------------------------------------------------------------------------------------------------------------------------- library(JuliaCall) library(tidyverse) library(tidygraph) @@ -10,664 +5,83 @@ library(ggraph) library(patchwork) library(pscl) -source("code/util.R") -source("code/expand_scenarios.R") - width <- 8 -height <- width/1.55 +height <- width/1.5 + +n_resample <- 250L JuliaCall::julia_setup(force = TRUE) JuliaCall::julia_source("code/setup.jl") +walk(list.files("code", pattern = "*.R", full.names = TRUE), source) -## ----plot-school-structure------------------------------------------------------------------------------------------------------------------------------------------------ - n_bubble <- 9 - bubbles_per_class <- 3 -classes_per_school <- 12 - - n_class <- n_bubble * bubbles_per_class -n_school <- classes_per_school * n_class - -pr_meet_class <- 3/(n_class - 1) -pr_meet_school <- 1/(n_school - 1) - -plt <- julia_call("school", - n_bubble, bubbles_per_class, classes_per_school, pr_meet_class, pr_meet_school, - 0.0, # gamma - 0.5, # frac symptomatic - 0.01, # pr_noncov_symptoms - Inf, # a for beta of test compliance - 1.0 # b - ) %>% - julia_call("get_adjacency_matrix", .) %>% - as_tbl_graph() %>% - ggraph("matrix") + - geom_edge_tile(aes(fill = weight)) + - scale_edge_fill_continuous( - "", - low = "#dddddd", high = "#000000", na.value = "white", - limits = c(0, NA_real_) - ) + - coord_fixed(expand = FALSE) + - theme( - legend.position = "right" - ) + - labs( - caption = "Adjacency matrix of school population;\r\nconnection strength is 'expected pairwise daily risk contacts' (across all shared groups)." - ) -save_plot(plt, "adjacency-matrix", width = 2*width/3, height = 2*height/3) - - +# calibrate-infectivity (standard case) ---------------------------------------- +source("plots/population-structure.R") -## ----calibrate-infectivity------------------------------------------------------------------------------------------------------------------------------------------------ -tbl_r_zero <- tibble( - gamma = seq(0, 0.1, length.out = 1000), - R = julia_call("simulate_r_zero.", - julia_call("school.", - n_bubble, bubbles_per_class, classes_per_school, pr_meet_class, pr_meet_school, - gamma, - 0.5, # frac symptomatic - 0.01, # pr_noncov_symptoms - Inf, 1.0 - ) - ) -) - -fit <- zeroinfl(formula = R ~ gamma | gamma, dist = "negbin", data = tbl_r_zero) - -# define inverse function -get_gamma <- function(R) { - f <- function(gamma) as.numeric(predict(fit, newdata = tibble(gamma = gamma), type = "response")) - as.numeric(uniroot(function(x) f(x) - R, interval = c(0, .1))$root) -} - -plt <- ggplot(tbl_r_zero) + - aes(gamma, R) + - geom_point(alpha = 0.2, shape = 16) + - geom_line( - data = tibble( - gamma = seq(0, .1, length.out = 100), - R = as.numeric(predict(fit, newdata = tibble(gamma = gamma), type = "response")) - ) - ) + - labs( - x = expression(gamma), - caption = "Fitted zero-inflated negative binomial regression for r-zero vs. infectivity." - ) -save_plot(plt, "calibrating-infectivity", width = width, height = height/2) +# calibrate-infectivity (standard case) ---------------------------------------- +source("plots/calibration-infectivity.R") -## ----plot-vl-trajectories, warning=FALSE---------------------------------------------------------------------------------------------------------------------------------- -dm <- julia_call("LarremoreModel", get_gamma(R = 3), frac_symptomatic = 0.5, need_return = "Julia") -individuals <- julia_call("Individual.", dm, rep(.01, 32), need_return = "Julia") -julia_call("infect!.", individuals, need_return = "Julia") -julia_call("steps!.", individuals, 21L, need_return = "Julia") +# plot viral load sample trajectories ------------------------------------------ +source("plots/vl-trajectories.R") -plt <- julia_call("get_status_logs", individuals, need_return = "R") %>% - as_tibble() %>% - mutate( - symptomatic = factor(symptomatic, levels = c(FALSE, TRUE), labels = c("non-symptomatic", "symptomatic")) - ) %>% - ggplot() + - aes(day, viral_load) + - geom_hline(yintercept = 1e6) + - geom_bar(aes(fill = symptomatic), stat = "identity", alpha = .66, width = 1) + - scale_y_log10() + - scale_x_continuous("day post infection", breaks = seq(0, 35, by = 7)) + - scale_fill_manual(breaks = c("non-symptomatic", "symptomatic"), values = c("darkgray", "black")) + - ylab("viral load") + - facet_wrap(~uuid, nrow = 4) + - theme( - strip.text = element_blank(), - strip.background = element_blank() - ) -save_plot(plt, "vl-trajectories") +# recalibrate innova lfd mean sensitivity -------------------------------------- +source("plots/fit-mean-sensitivity.R") -## ----define-pcr-test------------------------------------------------------------------------------------------------------------------------------------------------------ -pcr <- julia_call("FixedTest", "pcr", 0.975, 1.0, lod = 300.0, need_return = "Julia") +# plot autocorrelation --------------------------------------------------------- +source("plots/autocorrelation.R") -## ----fit-innova-data------------------------------------------------------------------------------------------------------------------------------------------------------ -# screen grabbed pixels -tbl_innova_data <- tibble( - `viral load` = 10^(2:7), - sensitivity = 1 - (c(621, 521, 344, 207, 144, 123) - 114) / (684 - 114) -) +# plot sensitivity vs infection probability ------------------------------------ +source("plots/sensitivity-vs-infectivity.R") -sensitivity <- function(vl, slope, intercept) { - lfd <- julia_call("LogRegTest", "lfd", slope, intercept, 0.998, need_return = "Julia") - julia_call("sensitivity.", lfd, vl, need_return = "R") -} - -innova <- optim( - c(.76, -5), - function(x) { - sum(( - tbl_innova_data$sensitivity - - sensitivity(tbl_innova_data$`viral load`, slope = x[1], intercept = x[2]) - )^2) - }, - method = "L-BFGS-B", - lower = c(.1, -10), - upper = c(2.5, 5), - control = list(maxit = 1e4) -) -innova$slope <- innova$par[1] -innova$intercept <- innova$par[2] - -## ----recalibrate-innova-lfd-test------------------------------------------------------------------------------------------------------------------------------------------ -dm <- julia_call("LarremoreModel", get_gamma(R = 3), frac_symptomatic = 0.5, need_return = "Julia") -individuals <- julia_call("Individual.", dm, rep(.01, 1e5), need_return = "Julia") -julia_call("infect!.", individuals, need_return = "Julia") -julia_call("steps!.", individuals, 21L, need_return = "Julia") - -tbl_vl <- julia_call("get_status_logs", individuals) %>% - as_tibble() %>% - arrange(uuid, day) %>% - group_by(uuid) %>% - filter( - row_number() >= which(viral_load > 300)[1], - row_number() < which(symptomatic)[1] - ) %>% - sample_n(1) %>% - ungroup() %>% - select(uuid, day, viral_load) - -mean_sensitivity <- function(eta) { - tbl_vl %>% - mutate( - sensitivity = sensitivity(viral_load^eta, slope = innova$slope, intercept = innova$intercept) - ) %>% - pull(sensitivity) %>% - summary() %>% - mean() -} - -# define inverse -get_eta <- function(sens) uniroot(function(x) mean_sensitivity(x) - sens, interval = c(0, 5))$root - -tbl_scenarios <- tibble( - `mean sensitivity` = c(mean_sensitivity(1), .4, .6, .8), - eta = map_dbl(`mean sensitivity`, get_eta) - ) - -plt <- tbl_scenarios %>% - expand_grid( - `viral load` = 10^seq(0, 11, length.out = 1000), - ) %>% - mutate( - sensitivity = sensitivity(`viral load`^eta, slope = innova$slope, intercept = innova$intercept) - ) %>% - ggplot() + - aes(`viral load`, sensitivity) + - geom_line(aes(linetype = sprintf("%.2f", `mean sensitivity`))) + - geom_point(data = tbl_innova_data) + - scale_x_log10() + - scale_y_continuous("", limits = c(0, 1)) + - scale_linetype_manual("", - breaks = sprintf("%.2f", tbl_scenarios$`mean sensitivity`), - values = c(1, 4, 3, 2), - labels = purrr::map2( - sprintf("%.2f", tbl_scenarios$`mean sensitivity`), - tbl_scenarios$eta, - ~bquote( .(..1) (eta == .(round(..2, 2))) ) - ) - ) + - labs( - caption = "Fitted Innova LFD test sensitivity." - ) + - theme( - legend.position = "right" - ) -save_plot(plt, "fit-sensitivity", width = width, height = height/2) - - - -## ----plot-acf------------------------------------------------------------------------------------------------------------------------------------------------------------- -get_ar <- function(ar_coefficient = 0, ar_window = 0L, seed = 42L, nrsmpl = 1e4) { - julia_call("Random.seed!", as.integer(seed)) - lfd <- julia_call("LogRegTest", - "lfd", innova$slope, innova$intercept, 0.998, ar_window = ar_window, ar_coefficient = ar_coefficient, - need_return = "Julia" - ) - dm <- julia_call("LarremoreModel", get_gamma(R = 3), frac_symptomatic = 0.5, need_return = "Julia") - individuals <- julia_call("Individual.", dm, rep(.01, nrsmpl), need_return = "Julia") - julia_call("infect!.", individuals, need_return = "Julia") - for (i in 1:21) { - if ((i %% 7) %in% 0:4) { # only test on weekends - julia_call("conduct_test!.", lfd, individuals, need_return = "None") - } - julia_call("step!.", individuals, need_return = "None") - } - tbl_status <- julia_call("get_status_logs", individuals, need_return = "R") %>% - as_tibble() - tbl_tests <- julia_call("get_test_logs", individuals, need_return = "R") %>% - as_tibble() %>% - select(-viral_load) - - left_join(tbl_tests, tbl_status, by = c("uuid", "day")) %>% - select(uuid, day, result, viral_load, probability_positive, symptomatic) %>% - group_by(uuid) %>% - filter( - !any(symptomatic) | (row_number() < which(symptomatic)[1]) - ) %>% - summarize( - acf = list(tibble( - lag = acf(as.numeric(result), plot = FALSE, lag.max = 7)$lag[,,1], - ac = acf(as.numeric(result), plot = FALSE, lag.max = 7)$acf[,,1] - )), - .groups = "drop" - ) %>% - unnest(acf) %>% - mutate( - ac = if_else(is.nan(ac), 1, ac) # constant signal has ac 1 - ) %>% - group_by(lag) %>% - summarize( - mid = mean(ac), - lo = quantile(ac, .05), - hi = quantile(ac, .95), - .groups = "drop" - ) -} - -plt <- tibble( - coefficient = c(0, .25, .5, .75) - ) %>% - mutate( - ar = map(coefficient, ~get_ar(., 3L, nrsmpl = 1e4)) - ) %>% - unnest(ar) %>% - mutate( - coefficient = sprintf("%.2f", coefficient) - ) %>% - ggplot() + - aes(lag, mid, color = coefficient) + - geom_line(aes(y = mid)) + - scale_color_grey("") + - ylab("autocorrelation") + - scale_x_continuous("lag [days]", breaks = seq(0, 21, by = 1)) + - theme( - legend.position = "right" - ) - -save_plot(plt, "test-autocorrelation", width = width, height = height/2) - - - -## ----plot-sensitivity-vs-infection-probability---------------------------------------------------------------------------------------------------------------------------- -plt <- expand_grid( - tibble( - `mean sensitivity` = c(.4, .6, .8), - eta = map_dbl(`mean sensitivity`, get_eta) +# define policies -------------------------------------------------------------- +pcr_turnaround <- 2L +isolation_duration <- 10L +lst_policies <- list( + reference = function(params) julia_call("SymptomaticIsolation", + lfd(params), pcr_test = pcr(params), + pcr_turnaround = pcr_turnaround, isolation_duration = isolation_duration, + need_return = "Julia" ), - tibble( - Rs = c(1.5, 3, 6), - gamma = map_dbl(Rs, get_gamma) - ) - ) %>% - expand_grid( - `viral load` = 10^seq(0, 11, length.out = 1000) - ) %>% - mutate( - sensitivity = sensitivity(`viral load`^eta, slope = innova$slope, intercept = innova$intercept), - `probability to infect` = pmax(0, pmin(1, gamma * (log10(`viral load`) - 6))) - ) %>% - ggplot() + - aes(`probability to infect`) + - geom_line(aes(y = sensitivity)) + - facet_grid( - `mean sensitivity` ~ Rs, - labeller = label_both - ) + - ylab("sensitivity") -save_plot(plt, "sensitivity-vs-infectivity", width = width, height = .9*height) - - - -## ----define-scenarios----------------------------------------------------------------------------------------------------------------------------------------------------- -tbl_sim_scenarios <- bind_rows( - # default scenario with different levels of symtoms - expand_scenarios( - frac_symptomatic = c(0.25, 0.5, 0.75) - ), - # different level of AR strength - expand_scenarios( - ar_coefficient = 0.75 - ), - # different bubble structure - expand_scenarios( - n_bubble = 27, - bubbles_per_class = 1 - ), - # different test compliance (mean = 66.67%) - expand_scenarios( - a = 2/15, b = 1/15 - ) -) - - -## ----run-simulation------------------------------------------------------------------------------------------------------------------------------------------------------- -tbl_results <- bind_cols( - tbl_sim_scenarios %>% - select(where(~length(unique(.)) > 1)), # drop constants - with(tbl_sim_scenarios, - julia_call("evaluate_performance", - n_bubble, bubbles_per_class, classes_per_school, pr_meet_class, pr_meet_school, - gamma, frac_symptomatic, pr_external_infections, pr_noncovid_symptoms, - pcr_test, - eta, slope, intercept, ar_window, ar_coefficient, - policy_name, - a, b, - need_return = "R" + `test for release` = function(params) julia_call("DynamicScreening", + lfd(params), pcr_test = pcr(params), + pcr_turnaround = pcr_turnaround, isolation_duration = isolation_duration, + followup_duration = 7L, + need_return = "Julia" + ), + `Thu/Fri off` = function(params) julia_call("SymptomaticIsolation", + lfd(params), pcr_test = pcr(params), + pcr_turnaround = pcr_turnaround, isolation_duration = isolation_duration, + fixed_isolation_weekdays = as.integer(c(3, 4)), + need_return = "Julia" + ), + `Mon screening` = function(params) julia_call("SymptomaticIsolation", + lfd(params), pcr_test = pcr(params), + pcr_turnaround = pcr_turnaround, isolation_duration = isolation_duration, + screening_test_weekdays = julia_eval("[0]", need_return = "Julia"), + need_return = "Julia" + ), + `Mon/Wed screening` = function(params) julia_call("SymptomaticIsolation", + lfd(params), pcr_test = pcr(params), + pcr_turnaround = pcr_turnaround, isolation_duration = isolation_duration, + screening_test_weekdays = as.integer(c(0, 2)), + need_return = "Julia" ) - ) -) %>% -rename( - `mean sensitivity` = mean_sensitivity -) %>% -mutate( - policy_name = str_replace(policy_name, "test & release", "test for release") ) -dir.create("_site") -readr::write_rds(tbl_results, "_site/tbl_results.rds", compress = "gz") - - -# main scenario ---------------------------------------------------------------- -plt <- tbl_results %>% - filter( - bubbles_per_class == 3, - ar_coefficient == 0.0, - frac_symptomatic == 0.5, - a == Inf, b == 1, - R == 3, - `mean sensitivity` == 0.6 - ) %>% - transmute( - policy_name, - `% infected (cumulative)` = n_infected/n_school, - `% schooldays missed (cumulative)` = workdays_missed/n_school/5/6 - ) %>% - pivot_longer(-policy_name) %>% - ggplot() + - aes(policy_name, value) + - geom_boxplot() + - scale_y_continuous("", labels = scales::percent) + - facet_wrap(~name) + - theme( - axis.text.x = element_text(angle = 33, hjust = 1), - axis.title.x = element_blank(), - legend.position = "right" - ) -save_plot(plt, "results-main", width = width, height = .75*height) - - - -## ------------------------------------------------------------------------------------------------------------------------------------------------------------------------- -plt <- tbl_results %>% - filter( - bubbles_per_class == 3, - ar_coefficient == 0.0, - frac_symptomatic == 0.5, - a == Inf, b == 1 - ) %>% - rename( - Rs = R - ) %>% - ggplot() + - aes(n_infected/n_school, mean_daily_infectious, color = policy_name, fill = policy_name) + - geom_vline(xintercept = 6/n_school) + - geom_point(alpha = .2, shape = 16) + - stat_ellipse(geom = "polygon", alpha = 0.33, color = NA, level = 0.9) + - scale_y_continuous("mean infectious (daily)") + - scale_x_continuous("% infected (cumulative)", labels = scales::percent, limits = c(NA, 1)) + - facet_grid(`mean sensitivity` ~ Rs, labeller = label_both) -save_plot(plt, "results-infectiousness-vs-infectivity", width = width, height = height) - - - -## ------------------------------------------------------------------------------------------------------------------------------------------------------------------------- -plt1 <- tbl_results %>% - filter( - bubbles_per_class == 3, - frac_symptomatic == 0.5, - a == Inf, b == 1 - ) %>% - rename( - Rs = R - ) %>% - ggplot() + - aes(policy_name, n_infected/n_school, color = factor(ar_coefficient)) + - geom_hline(yintercept = 6/n_school) + - geom_boxplot() + - scale_y_continuous("% infected (cumulative)", labels = scales::percent, limits = c(0, 1)) + - # scale_color_grey() + - facet_grid(`mean sensitivity` ~ Rs, labeller = label_both) + - theme( - axis.text.x = element_text(angle = 33, hjust = 1), - axis.title.x = element_blank(), - legend.position = "bottom" - ) + - ggtitle("Distribution of % infected by AR coefficient") -plt2 <- tbl_results %>% - filter( - bubbles_per_class == 3, - a == Inf, b == 1 - ) %>% - rename( - Rs = R - ) %>% - ggplot() + - aes(policy_name, n_infected/n_school, color = factor(1 - frac_symptomatic)) + - geom_hline(yintercept = 6/n_school) + - geom_boxplot() + - scale_y_continuous("% infected (cumulative)", labels = scales::percent, limits = c(0, 1)) + - # scale_color_grey() + - facet_grid(`mean sensitivity` ~ Rs, labeller = label_both) + - theme( - axis.text.x = element_text(angle = 33, hjust = 1), - axis.title.x = element_blank(), - legend.position = "bottom" - ) + - ggtitle("Distribution of % infected by % asymptomatic") -plt <- plt2 + plt1 + plot_layout(nrow = 1) -save_plot(plt, "results-infectivity-marginal", width = 1.33*width, height = 1.33*height) - - - -## ------------------------------------------------------------------------------------------------------------------------------------------------------------------------- -plt <- tbl_results %>% - filter( - bubbles_per_class == 3, - ar_coefficient == 0.0, - frac_symptomatic == 0.5, - a == Inf, b == 1 - ) %>% - rename( - Rs = R - ) %>% - ggplot() + - aes(policy_name, n_pcr_tests) + - geom_hline(yintercept = 6/n_school) + - geom_boxplot() + - scale_y_continuous("PCR test required (cumulative)") + - facet_grid(`mean sensitivity` ~ Rs, labeller = label_both) + - theme( - axis.text.x = element_text(angle = 33, hjust = 1), - axis.title.x = element_blank(), - legend.position = "right" - ) -save_plot(plt, "results-pcr-marginal", width = width, height = height) +# main results, base scenario -------------------------------------------------- +source("plots/main-results.R") +# only one bubble per class ---------------------------------------------------- +source("plots/only-one-bubble.R") +# less compliance -------------------------------------------------------------- +source("plots/lower-lfd-compliance.R") -## ------------------------------------------------------------------------------------------------------------------------------------------------------------------------- -plt <- tbl_results %>% - filter( - bubbles_per_class == 3, - ar_coefficient == 0.0, - frac_symptomatic == 0.5, - a == Inf, b == 1 - ) %>% - rename( - Rs = R - ) %>% - ggplot() + - aes(n_infected/n_school, n_pcr_tests, color = policy_name, fill = policy_name) + - geom_vline(xintercept = 6/n_school) + - geom_point(alpha = .2, shape = 16) + - stat_ellipse(geom = "polygon", alpha = 0.33, color = NA, level = 0.9) + - scale_y_continuous("PCR test required (cumulative)") + - scale_x_continuous("% infected (cumulative)", labels = scales::percent, limits = c(NA, 1)) + - facet_grid(`mean sensitivity` ~ Rs, labeller = label_both) -save_plot(plt, "results-pcr-vs-infectivity", width = width, height = height) - - - -## ------------------------------------------------------------------------------------------------------------------------------------------------------------------------- -plt <- tbl_results %>% - filter( - bubbles_per_class == 3, - ar_coefficient == 0.0, - frac_symptomatic == 0.5, - a == Inf, b == 1 - ) %>% - rename( - Rs = R - ) %>% - ggplot() + - aes(n_infected/n_school, workdays_missed/n_school/5/6, color = policy_name, fill = policy_name) + - geom_abline(slope = 1) + - geom_point(alpha = .2, shape = 16) + - stat_ellipse(geom = "polygon", alpha = 0.33, color = NA, level = 0.9) + - scale_x_continuous("% infected (cumulative)", labels = scales::percent, limits = c(NA, 1)) + - scale_y_continuous(labels = scales::percent) + - facet_grid(`mean sensitivity` ~ Rs, labeller = label_both) + - labs( - y = "% schooldays missed (cumulative)" - ) -save_plot(plt, "results-schooldays-missed-vs-infectivity", width = width, height = height) - - - -## ------------------------------------------------------------------------------------------------------------------------------------------------------------------------- -plt <- tbl_results %>% - filter( - bubbles_per_class == 3, - ar_coefficient == 0.75, - frac_symptomatic == 0.5, - a == Inf, b == 1 - ) %>% - rename( - Rs = R - ) %>% - ggplot() + - aes(n_infected/n_school, workdays_missed/n_school/5/6, color = policy_name, fill = policy_name) + - geom_abline(slope = 1) + - geom_point(alpha = .2, shape = 16) + - stat_ellipse(geom = "polygon", alpha = 0.33, color = NA, level = 0.9) + - scale_x_continuous("% infected (cumulative)", labels = scales::percent, limits = c(NA, 1)) + - scale_y_continuous(labels = scales::percent) + - facet_grid(`mean sensitivity` ~ Rs, labeller = label_both) + - labs( - y = "% schooldays missed (cumulative)" - ) -save_plot(plt, "sensitivity-schooldays-missed-vs-infectivity-ar", width = width, height = height) - - - -## ------------------------------------------------------------------------------------------------------------------------------------------------------------------------- -plt <- tbl_results %>% - filter( - bubbles_per_class == 1, - ar_coefficient == 0.0, - frac_symptomatic == 0.5, - a == Inf, b == 1 - ) %>% - rename( - `Rs*` = R - ) %>% - ggplot() + - aes(n_infected/n_school, workdays_missed/n_school/5/6, color = policy_name, fill = policy_name) + - geom_abline(slope = 1) + - geom_point(alpha = .2, shape = 16) + - stat_ellipse(geom = "polygon", alpha = 0.33, color = NA, level = 0.9) + - scale_x_continuous("% infected (cumulative)", labels = scales::percent, limits = c(NA, 1)) + - scale_y_continuous(labels = scales::percent) + - facet_grid(`mean sensitivity` ~ `Rs*`, labeller = label_both) + - labs( - y = "% schooldays missed (cumulative)" - ) -save_plot(plt, "sensitivity-schooldays-missed-vs-infectivity-bubbles", width = width, height = height) - -plt <- tbl_results %>% - filter( - ar_coefficient == 0.0, - frac_symptomatic == 0.5, - a == Inf, b == 1 - ) %>% - rename( - `Rs*` = R - ) %>% - ggplot() + - aes(policy_name, n_infected/n_school, color = factor(bubbles_per_class)) + - geom_hline(yintercept = 6/n_school) + - geom_boxplot() + - scale_y_continuous("% infected (cumulative)", labels = scales::percent, limits = c(0, 1)) + - scale_color_discrete("number of bubbles per class:") + - facet_grid(`mean sensitivity` ~ `Rs*`, labeller = label_both) + - theme( - axis.text.x = element_text(angle = 33, hjust = 1), - axis.title.x = element_blank(), - legend.title = element_text() - ) -save_plot(plt, "sensitivity-infectivity-marignal-bubbles", width = width, height = 1.1*height) - - -# sensitivity wrt to compliance ------------------------------------------------ -# plot beta distribution used to sample individual compliance -plt <- tibble( - compliance = rbeta(10^5, shape1 = 2/15, shape2 = 1/15) - ) %>% - ggplot() + - aes(compliance) + - geom_histogram(aes(y = ..ncount..), binwidth = 0.05) + - scale_y_continuous("", labels = scales::percent) -save_plot(plt, "sensitivity-infectivity-compliance-beta", width = .66*width, height = .5*height) - -# plot results -plt <- tbl_results %>% - filter( - bubbles_per_class == 3, - ar_coefficient == 0.0, - frac_symptomatic == 0.5 - ) %>% - rename( - Rs = R - ) %>% - mutate( - tmp = if_else(is.finite(a), a / (a + b), 1), - test_compliance = factor( - tmp, - levels = unique(tmp), - labels = sprintf("%.2f", unique(tmp)) - ) - ) %>% - ggplot() + - aes(policy_name, n_infected/n_school, color = test_compliance) + - geom_hline(yintercept = 6/n_school) + - geom_boxplot() + - scale_y_continuous("% infected (cumulative)", labels = scales::percent, limits = c(0, 1)) + - scale_color_discrete("average LFD test compliance probability:") + - facet_grid(`mean sensitivity` ~ Rs, labeller = label_both) + - theme( - axis.text.x = element_text(angle = 33, hjust = 1), - axis.title.x = element_blank(), - legend.title = element_text() - ) -save_plot(plt, "sensitivity-infectivity-compliance", width = width, height = 1.1*height) +# lower limit of infectivity --------------------------------------------------- +source("plots/lower-lli.R")