From aeda465208168fac67b7a497ac60a1d295688afd Mon Sep 17 00:00:00 2001 From: "Carl A. B. Pearson" Date: Fri, 13 Jan 2023 07:21:26 +0000 Subject: [PATCH] visualization updates --- exp/active-vac/fig/Makefile | 8 +-- exp/active-vac/fig/alt_i_all.R | 64 ++++++++++++++++++++++++ exp/active-vac/fig/alt_inc_non.R | 86 ++++++++++++++++++++++++++++++++ exp/active-vac/fig/vis_support.R | 12 ++--- 4 files changed, 160 insertions(+), 10 deletions(-) create mode 100644 exp/active-vac/fig/alt_i_all.R create mode 100644 exp/active-vac/fig/alt_inc_non.R diff --git a/exp/active-vac/fig/Makefile b/exp/active-vac/fig/Makefile index bbf38ae..ee7a5db 100644 --- a/exp/active-vac/fig/Makefile +++ b/exp/active-vac/fig/Makefile @@ -82,10 +82,6 @@ ${OUTDIR}/cum_averted_%.png: cum_avert_outcome.R vis_support.rda ${PLOTINS} | ${ ${OUTDIR}/cum_eff_%.png: cum_eff_outcome.R vis_support.rda ${PLOTINS} | ${OUTDIR} $(call R) -${OUTDIR}/alt_eff_%.png: alt_eff_outcome.R vis_support.rda ${RESDIR}/alt_eff.rds \ -${RESDIR}/digest-key.rds | ${OUTDIR} - $(call R) - ${OUTDIR}/alt_ave_%.png: alt_ave_outcome.R vis_support.rda ${RESDIR}/alt_eff.rds \ ${RESDIR}/digest-key.rds | ${OUTDIR} $(call R) @@ -110,6 +106,10 @@ ${OUTDIR}/alt_inc_%.png: alt_inc_outcome.R vis_support.rda ${RESDIR}/alt_eff.rds ${RESDIR}/digest-key.rds | ${OUTDIR} $(call R) +${OUTDIR}/alt_eff_%.png: alt_eff_outcome.R vis_support.rda ${RESDIR}/alt_eff.rds \ +${RESDIR}/digest-key.rds | ${OUTDIR} + $(call R) + allavert: $(addprefix ${OUTDIR}/,$(patsubst %,avertvis_%.png,inf sev deaths) $(patsubst %,cum_averted_%.png,inf sev deaths) avertvis.png) alleff: $(addprefix ${OUTDIR}/,$(patsubst %,cum_eff_%.png,inf sev deaths)) alteff: $(addprefix ${OUTDIR}/,$(patsubst %,alt_eff_%.png,inf sev deaths)) diff --git a/exp/active-vac/fig/alt_i_all.R b/exp/active-vac/fig/alt_i_all.R new file mode 100644 index 0000000..801a557 --- /dev/null +++ b/exp/active-vac/fig/alt_i_all.R @@ -0,0 +1,64 @@ + +.pkgs <- c("data.table", "ggplot2", "scales", "ggh4x", "cabputils", "geomtextpath") +.pkgs |> sapply(require, character.only = TRUE) |> all() |> stopifnot() + +#' assumes R project at the experiment root level +.args <- if (interactive()) c( + file.path("fig", "vis_support.rda"), + file.path("fig", "process", c("digest.rds", "digest-key.rds", "vocwindows.rds")), + file.path("fig", "output", "alt_inc_all.png") +) else commandArgs(trailingOnly = TRUE) + +load(.args[1]) + +intfilter <- expression(realization >= 0) + +#' comes key'd +inc.dt <- readRDS(.args[2])[ + eval(datefilter) & eval(intfilter) & eval(outfilter) +][, .( + scenario, realization, date, outcome, value +)] + +intscns <- inc.dt[, unique(scenario)] + +scn.dt <- readRDS(.args[3]) + +takeover.wins <- readRDS(.args[4]) + +intscn.dt <- scn.dt[scenario %in% intscns] +# reconstructing reference scenarios +refscn.dt <- scn.dt[quar == FALSE & pas_vac == TRUE & act_vac == "none"] + +# incref.dt <- inc.dt[ +# intscn.dt, on=.(scenario) +# ][(act_vac == "ring") & (quar == FALSE)][ # only need to go from one reference +# refscn.dt, on =.(act_alloc = pas_alloc, inf_con) +# ][,.( +# value = value[1] +# ), by=.(scenario = i.scenario, realization, date, outcome) +# ][refscn.dt, on=.(scenario)] + +plt.dt <- setkeyv( + inc.dt[intscn.dt, on=.(scenario)], + union(key(inc.dt), colnames(scn.dt)) +)[inf_con == FALSE] + +rm(inc.dt) +gc() + +plt.qs <- plt.prep(plt.dt, j = .(value)) + +mm.ref <- plt.qs[,.(ymin = min(q90l), ymax = max(q90h)),by=.(outcome)] +mm.ref[, ymin := ymin - .15*(ymax-ymin) ] +mm.ref[, yspan := ymax - ymin ] +tw <- takeover.wins[q == 0.5][CJ(measure, outcome = mm.ref$outcome), on=.(measure)][mm.ref, on=.(outcome)] +tw[, end := pmin(end, plt.qs[, max(date)])] +tw[, mid := start + (end-start)/2 ] + +p <- allplot( + plt.qs, yl = "Per 10k, Incidence of ...", + withRef = FALSE, ins = voc.band(tw) +) + +ggsave(tail(.args, 1), p, height = 6, width = 10, bg = "white") diff --git a/exp/active-vac/fig/alt_inc_non.R b/exp/active-vac/fig/alt_inc_non.R new file mode 100644 index 0000000..0a41e42 --- /dev/null +++ b/exp/active-vac/fig/alt_inc_non.R @@ -0,0 +1,86 @@ + +.pkgs <- c("data.table", "ggplot2", "scales", "ggh4x", "cabputils", "geomtextpath") +.pkgs |> sapply(require, character.only = TRUE) |> all() |> stopifnot() + +#' assumes R project at the experiment root level +.args <- if (interactive()) c( + file.path("fig", "vis_support.rda"), + file.path("fig", "process", c("digest-ref.rds", "digest-key.rds", "vocwindows.rds", "digest.rds")), + file.path("fig", "output", "alt_inc_non.png") +) else commandArgs(trailingOnly = TRUE) + +load(.args[1]) + +intfilter <- expression(realization >= 0) + +#' comes key'd - all the intervention results +inc.dt <- readRDS(.args[2])[ + eval(datefilter) & eval(intfilter) & eval(outfilter) +][, .( + scenario, realization, date, outcome, value +)] + +scn.dt <- readRDS(.args[3]) +takeover.wins <- readRDS(.args[4]) +# reconstructing reference scenarios +refscn.dt <- scn.dt[(quar == FALSE & act_vac == "none") & (!inf_con | !pas_vac)] + +int.dt <- readRDS(.args[5])[ + eval(datefilter) & eval(intfilter) & eval(outfilter) +][scenario %in% refscn.dt$scenario, .( + scenario, realization, date, outcome, value +)] + +plt.dt <- setkeyv( + rbind(inc.dt, int.dt)[refscn.dt, on=.(scenario)], + union(key(inc.dt), colnames(scn.dt)) +) + +rm(inc.dt) +gc() + +plt.qs <- plt.prep(plt.dt, j = .(value)) + +mm.ref <- plt.qs[,.(ymin = min(q90l), ymax = max(q90h)),by=.(outcome)] +mm.ref[, ymin := ymin - .15*(ymax-ymin) ] +mm.ref[, yspan := ymax - ymin ] +tw <- takeover.wins[q == 0.5][CJ(measure, outcome = mm.ref$outcome), on=.(measure)][mm.ref, on=.(outcome)] +tw[, end := pmin(end, plt.qs[, max(date)])] +tw[, mid := start + (end-start)/2 ] + +allplot <- function (data.qs, yl, withRef = FALSE, col.breaks = if (withRef) c("risk", + "age", "ring") else c("none", "risk", "age", "ring"), withBands = NULL, + ins = list()) +{ + res <- ggplot(data.qs) + aes( + x = date, color = pas_alloc, linetype = factor(c("nonpi", + "wquar")[quar + 1])) + facet_nested(rows = vars(outcome), + switch = "y", scales = "free_y", + labeller = labeller(outcome = c(inf = "Infection", sev = "Severe Disease", + deaths = "Deaths", vaccine = "Per 10K,\nDoses Administered"))) + + geom_month_background(data.qs, by = c(row = "outcome", + col = "talloc"), font.size = 3, value.col = "qmed", + max.col = "q90h", min.col = "q90l") + ins + if (withRef) { + res <- res + geom_hline(aes(yintercept = 0, color = "none", + linetype = "nonpi"), show.legend = FALSE, data = function(dt) dt[, + .SD[1], by = .(outcome, talloc)]) + } + res <- res + geom_ribbon(aes(ymin = q90l, ymax = q90h, fill = pas_alloc, + color = NULL), alpha = 0.15) + geom_line(aes(y = qmed)) + + scale_color_discrete(name = "Supply", aesthetics = c("color", "fill")) + + coord_cartesian(clip = "off") + + scale_y_continuous(name = yl, expand = c(0, 0), labels = scales::label_number(scale_cut = scales::cut_short_scale())) + + scale_x_null() + scale_linetype_quar(guide = "none") + + theme_minimal() + + theme(legend.position = "bottom", strip.placement = "outside", + legend.direction = "horizontal") + return(res) +} + +p <- allplot( + plt.qs, yl = "Per 10k, Incidence of ...", + ins = voc.band(tw) +) + +ggsave(tail(.args, 1), p, height = 6, width = 10, bg = "white") diff --git a/exp/active-vac/fig/vis_support.R b/exp/active-vac/fig/vis_support.R index 8a3be7a..1c897f4 100644 --- a/exp/active-vac/fig/vis_support.R +++ b/exp/active-vac/fig/vis_support.R @@ -561,12 +561,12 @@ allplot <- function( geom_hline( aes(yintercept = 0, color = "none", linetype = "nonpi"), show.legend = FALSE, data = \(dt) dt[,.SD[1],by=.(outcome, talloc)] - ) + - geom_texthline( - aes(yintercept = 0, color = "none", label = "Reference\nProgram"), - inherit.aes = FALSE, show.legend = FALSE, data = \(dt) dt[talloc == "LIC",.SD[1],by=.(outcome, talloc)], - hjust = 0, gap = FALSE - ) + )# + + # geom_texthline( + # aes(yintercept = 0, color = "none", label = "Reference\nProgram"), + # inherit.aes = FALSE, show.legend = FALSE, data = \(dt) dt[talloc == "LIC",.SD[1],by=.(outcome, talloc)], + # hjust = 0, gap = FALSE + # ) } res <- res + geom_ribbon(aes(ymin=q90l, ymax=q90h, fill=act_vac, color=NULL), alpha=0.15) +