From 7ec5f62ec5420eb50a7a89b850920abb45acdcd9 Mon Sep 17 00:00:00 2001 From: "Carl A. B. Pearson" Date: Fri, 31 Jul 2020 08:24:49 -0400 Subject: [PATCH] detailed PAK scenarios --- setup/PAK_detail/Makefile | 41 ++++ setup/PAK_detail/create_params_set.R | 186 ++++++++++++++++++ .../PAK_detail/generate_scenario_parameters.R | 103 ++++++++++ simulate/run_scenarios_detail.R | 27 +-- 4 files changed, 337 insertions(+), 20 deletions(-) create mode 100644 setup/PAK_detail/Makefile create mode 100644 setup/PAK_detail/create_params_set.R create mode 100644 setup/PAK_detail/generate_scenario_parameters.R diff --git a/setup/PAK_detail/Makefile b/setup/PAK_detail/Makefile new file mode 100644 index 0000000..8014f3e --- /dev/null +++ b/setup/PAK_detail/Makefile @@ -0,0 +1,41 @@ + +-include ../../local.makefile + +COVIDMPATH ?= ../../../covidm +DATADIR ?= ~/Dropbox/covidm_reports +#DATADIR ?= ~/Dropbox/covidm_test +SAMPN ?= 500 +INDIR := ${DATADIR}/generation +OUTDIR := ${DATADIR}/hpc_inputs +DETDIR := ${DATADIR}/hpc_detailed + +R = Rscript $^ $@ + +# TODO other directory structure dependencies +${DETDIR}: + mkdir -p $@ + +params: create_params_set.R ${DETDIR}/worldpop5yr.lfs.csv + Rscript $^ ${COVIDMPATH} ${DETDIR} + +# currently, no explicit province contact matrices, so use national one +${DETDIR}/contact_matrices.rds: ${OUTDIR}/PAK/contact_matrices.rds + cp $< $@ + +# TODO update dependencies +# ${DETDIR}/%/params_set.rds: create_params_set.R ${INDIR}/data_contacts_missing.csv +# mkdir -p $(@D) +# Rscript $^ ${COVIDMPATH} $* $@ + +${DETDIR}/timing.rds: ${OUTDIR}/PAK/timing.rds + cp $< $@ + +${DETDIR}/run_options.rds: ${OUTDIR}/run_options.rds + cp $< $@ + +${DETDIR}/alt_scenarios.rds: generate_scenario_parameters.R + Rscript $^ $@ + +${DETDIR}/report_ref.rds: ${DETDIR}/alt_scenarios.rds + +dosetup: params ${DETDIR}/alt_scenarios.rds \ No newline at end of file diff --git a/setup/PAK_detail/create_params_set.R b/setup/PAK_detail/create_params_set.R new file mode 100644 index 0000000..1c5031b --- /dev/null +++ b/setup/PAK_detail/create_params_set.R @@ -0,0 +1,186 @@ +#' default simulation parameters +suppressPackageStartupMessages({ + require(data.table) +}) + +.args <- if (interactive()) c( + "~/Dropbox/covidm_reports/hpc_detailed/worldpop5yr.lfs.csv", + "../covidm", + "~/Dropbox/covidm_reports/hpc_detailed" # going to be writing admin#/param_set.rds +) else commandArgs(trailingOnly = TRUE) +#' @examples +#' .args <- gsub("ZWE","guineabissau",.args) +#' .args <- gsub("ZWE","palestine",.args) + +reference = fread( + .args[1], strip.white = FALSE)[ + level == 1 & country == "PAK" +] +country <- matref <- "Pakistan" + +cm_path = .args[2] +outfilefmt = file.path(tail(.args, 1),"%s","params_set.rds") + +cm_force_rebuild = F; +cm_build_verbose = F; +cm_force_shared = T; +cm_version = 2; + +suppressPackageStartupMessages({ + source(file.path(cm_path, "R", "covidm.R")) +}) + +# country <- cm_populations[ +# country_code == popcode, +# unique(as.character(name)) +# ] + +#set up node-runs +probs = fread( + "Age_low,Age_high,Prop_symptomatic,IFR,Prop_inf_hosp,Prop_inf_critical,Prop_critical_fatal,Prop_noncritical_fatal,Prop_symp_hospitalised,Prop_hospitalised_critical +0,9,0.66,8.59E-05,0.002361009,6.44E-05,0.5,0,0,0.3 +10,19,0.66,0.000122561,0.003370421,9.19E-05,0.5,9.47E-04,0.007615301,0.3 +20,29,0.66,0.000382331,0.010514103,0.000286748,0.5,0.001005803,0.008086654,0.3 +30,39,0.66,0.000851765,0.023423527,0.000638823,0.5,0.001231579,0.009901895,0.3 +40,49,0.66,0.001489873,0.0394717,0.001117404,0.5,0.002305449,0.018535807,0.3 +50,59,0.66,0.006933589,0.098113786,0.005200192,0.5,0.006754596,0.054306954,0.3 +60,69,0.66,0.022120421,0.224965092,0.016590316,0.5,0.018720727,0.150514645,0.3 +70,79,0.66,0.059223786,0.362002579,0.04441784,0.5,0.041408882,0.332927412,0.3 +80,100,0.66,0.087585558,0.437927788,0.065689168,0.5,0.076818182,0.617618182,0.3" +) + +#increase CFR +cfr_RR <- 1.5 +probs[, Prop_critical_fatal := Prop_critical_fatal * cfr_RR] +probs[, Prop_noncritical_fatal := Prop_noncritical_fatal * cfr_RR] + +#min(1, x) does not work for noncritical_fatal, for some reason +probs[Prop_critical_fatal > 1, Prop_critical_fatal := 1] +probs[Prop_noncritical_fatal > 1, Prop_noncritical_fatal := 1] + +reformat = function(P, lmic_adjust=TRUE) { + # no info to re-weight these, so assume 70-74 is like 70-79, and 75+ is like 80+ + if(lmic_adjust){ + P <- P[2:length(P)] + return (rep(P[1:8], each = 2)) + } else { + return (c(rep(P[1:7], each = 2),P[8:9])) + } +} + +P.icu_symp = reformat(probs[, Prop_symp_hospitalised * Prop_hospitalised_critical], lmic_adjust = TRUE); +P.nonicu_symp = reformat(probs[, Prop_symp_hospitalised * (1 - Prop_hospitalised_critical)], lmic_adjust = TRUE); +P.death = reformat(probs[, Prop_noncritical_fatal], lmic_adjust = TRUE); + +max_time <- 60 +tres <- 0.25 + +ponset2hosp <- cm_delay_gamma(7, 7, max_time, tres)$p +pignore <- cm_delay_skip(max_time, tres)$p +icustay <- cm_delay_gamma(10, 10, max_time, tres)$p +nonicustay <- cm_delay_gamma(8, 8, max_time, tres)$p +ponset2death <- cm_delay_gamma(22, 22, max_time, tres)$p + +cm_multinom_process <- function( + src, outcomes, delays, + report = "" +) { + if ("null" %in% names(outcomes)) { + if (length(report) != length(outcomes)) report <- rep(report, length(outcomes)) + report[which(names(outcomes)=="null")] <- "" + if (!("null" %in% names(delays))) { + delays$null <- c(1, rep(0, length(delays[[1]])-1)) + } + } else if (!all(rowSums(outcomes)==1)) { + report <- c(rep(report, length(outcomes)), "") + outcomes$null <- 1-rowSums(outcomes) + delays$null <- c(1, rep(0, length(delays[[1]])-1)) + } + nrow <- length(outcomes) + list( + source = src, type="multinomial", names=names(outcomes), report = report, + prob = t(as.matrix(outcomes)), delays = t(as.matrix(delays)) + ) +} + +cm_track_process <- function(src, name, delays, agecats = 16, report = "p") { + list( + source = src, type="multinomial", names = name, report = report, + prob = matrix(1, nrow = 1, ncol = agecats), + delays = t(delays) + ) +} + +burden_processes = list( + # process of sending symptomatic cases to hospital icu or ward + cm_multinom_process( + "Ip", + outcomes = data.frame(to_icu = P.icu_symp, to_nonicu = P.nonicu_symp), + delays = data.frame(to_icu = ponset2hosp, to_nonicu = ponset2hosp) + ), + # track icu prevalance + cm_track_process("to_icu", "icu", icustay), + # track ward prevalence + cm_track_process("to_nonicu", "nonicu", nonicustay), + # track infections - get from delta R prevalence +# cm_track_process("S", "infection", pignore, report="i"), + # send some cases to death, tracking outcidence + cm_multinom_process( + "Ip", + outcomes = data.table(death=P.death), + delays = data.table(death=ponset2death), + report = "o" + ) +) + +popnorm <- function(x, seed_cases = 50){ + #no contact matrix - will be overwritten by empirical ones + x$matrices$home <- x$matrices$work <- x$matrices$school <- x$matrices$other <- x$matrices$other*0 + + #age-specific probability of being symptomatic + #x$y <- c(rep(0.056, 3), rep(0.49, 8), rep(0.74, 8)) + #new values proposed by nick + x$y <- c( + rep(0.2973718, 2), rep(0.2230287, 2), rep(0.4191036, 2), + rep(0.4445867, 2), rep(0.5635720, 2), rep(0.8169443, 6) + ) + + #no cases in empty compartments + x$dist_seed_ages <- as.numeric(!(x$size == 0)) + + #seed cases + x$seed_times <- rep(0, seed_cases) + + return(x) +} + +for (admin_code in reference[, key]) { + params_set <- list() + + params1 <- cm_parameters_SEI3R( + country, matref, + deterministic=FALSE, date_end = 365*2 + ) + + extractpop <- melt(reference[ + key == admin_code, .SD, .SDcols = grep("(f|m)_", colnames(reference), value = TRUE) + ], measure.vars = grep("(f|m)_", colnames(reference), value = TRUE)) + extractpop[, agelb := as.integer(gsub(".+_(\\d+)","\\1", variable)) ] + redage <- extractpop[,.(value = sum(value)), keyby=agelb ] + redage[16]$value <- sum(redage[16:.N]$value) + + params1$pop[[1]]$size <- redage[1:16, value] + + params1$processes = burden_processes + + params1$pop <- lapply(params1$pop, popnorm) + #params1$time1 <- as.Date(params1$time1) + + params_set[[1]] <- params1 + + dir.create(dirname(sprintf(outfilefmt, admin_code))) + saveRDS(params_set, sprintf(outfilefmt, admin_code)) + +} + +#general population parameters diff --git a/setup/PAK_detail/generate_scenario_parameters.R b/setup/PAK_detail/generate_scenario_parameters.R new file mode 100644 index 0000000..521b3a1 --- /dev/null +++ b/setup/PAK_detail/generate_scenario_parameters.R @@ -0,0 +1,103 @@ +suppressPackageStartupMessages({ + require(data.table) +}) + +.args <- if (interactive()) c( + "~/Dropbox/covidm_reports/hpc_detailed/alt_scenarios.rds" +) else commandArgs(trailingOnly = TRUE) + +outfile <- tail(.args, 1) +reffile <- gsub("/(\\w+)\\.rds", "/report_ref.rds", outfile) + +# say that 11 March was the 50 active infections date +mar23 <- as.integer(as.Date("2020-03-23") - as.Date("2020-03-11")) + +dg <- function(...) data.table(expand.grid(..., stringsAsFactors = FALSE)) + +ref <- function( + age_split = NA, + self_iso = 0, + population = -1, coverage = 1, # coverage generally only applies to multi-pop models + school = 0, home = 0, work = 0, other = 0, # recall == these are reductions + travel = 0, + start_day = NA_integer_, end_day = Inf, + trigger_type = NA_character_, + trigger_value = NA_real_, + ... +) do.call(dg, c(as.list(environment()), list(...))) + +#' TODO: ignore? going to just iterate over rows where scen_id = 1, and if it's empty, nothing to do +unmitigated <- ref(scen_type = "unmitigated")[, scen_id := 1 ] + +reportref <- data.table(scen_id=integer(), label=character()) +append_scenarios <- function(dt, new_scen_ids, new_labels, scen_type) { + newdt <- data.table(scen_id = new_scen_ids, label = new_labels) + rbind(dt, newdt) +} + +scen_counter <- 1 +tagscenario <- function(dt, sc) { + dt[, scen_id := (1:.N) + sc ] + return(dt[.N, scen_id]) +} + +# cat(sprintf("unmitigated: %i", 1)) +reportref <- append_scenarios(reportref, 1L, "Unmitigated") + +#' all remaining scenarios consider some level of self-isolation with symptoms +#' and potentially different school effects +self_iso_lvls <- c(0.25) +schoolslowhi <- c(low=0.8, hi=1) # on = 0 reduction, off = 100% reduction + +#' - generic social distancing in work / other, levels 20-80%, by 10% +#' social distancing levels; use minimum for combination interventions + +work_dist <- c(0.2, 0.4) +home_dist <- c(0.2, 0.4) +other_dist <- c(0, 0.2) + +low_soc_all <- ref( + self_iso = self_iso_lvls[1], + home = home_dist, + school = schoolslowhi, + other = other_dist, + work = work_dist, + trigger_type = "day", + trigger_value = mar23 +)[, start_day := 0 ][, end_day := 90 ] +scen_counter <- tagscenario(low_soc_all, scen_counter) + +# cat(sprintf("40%% elder sheilding at home + low level all other measures: %i", elder_shielding[home == 0.4, scen_id])) +reportref <- append_scenarios( + reportref, + low_soc_all[, scen_id], + sprintf("generic options") +) + +#' - schools shutdown + continuous low level social distance (20%) + weekly alternating extra social distancing (40-50-60-70) => that + 20% +alternating <- ref( + self_iso = self_iso_lvls[1], + school = schoolslowhi["hi"], + home = max(home_dist), + work = max(work_dist), + other = max(other_dist), + start_day = mar23, end_day = mar23+30, + trigger_type = "stride", + trigger_value = 90 +) +scen_counter <- tagscenario(alternating, scen_counter) + +# cat("no school + intermittent 40%% lockdown: %i", alternating[work == 0.5, scen_id]) +reportref <- append_scenarios( + reportref, + alternating[, scen_id], + sprintf("alternating") +) + +all_scenarios <- setkey(rbind( + low_soc_all, + alternating +), scen_id) + +saveRDS(reportref, reffile) +saveRDS(all_scenarios, outfile) diff --git a/simulate/run_scenarios_detail.R b/simulate/run_scenarios_detail.R index 2567d26..5c8e7bc 100644 --- a/simulate/run_scenarios_detail.R +++ b/simulate/run_scenarios_detail.R @@ -4,9 +4,9 @@ suppressPackageStartupMessages({ }) .args <- if (interactive()) c( - "../covidm", "ZWE", "001", - sprintf("~/Dropbox/covidm_reports/hpc_inputs"), - "simulate/ZWE/001.qs" + "../covidm", "2428", "001", + sprintf("~/Dropbox/covidm_reports/hpc_detailed"), + "simulate/detail/001.qs" ) else commandArgs(trailingOnly = TRUE) # load covidm @@ -14,7 +14,7 @@ cm_path = .args[1] cm_force_rebuild = F; cm_build_verbose = F; cm_force_shared = T -cm_version = 1 +cm_version = 2 suppressPackageStartupMessages({ source(paste0(cm_path, "/R/covidm.R")) @@ -102,17 +102,6 @@ uref <- unname(as.matrix(yu[, .SD, .SDcols = grep("u_",colnames(yu))])) has_age_split <- iv_data[, any(!is.na(age_split))] -reduce_ages <- function (dt) { - fctr <- function(i, lvls = c("<14", "15-29", "30-44","45-59", "60+")) factor( - lvls[i], levels = lvls, ordered = T - ) - dt[between(as.integer(group), 1, 3), age := fctr(1) ] - dt[between(as.integer(group), 4, 6), age := fctr(2) ] - dt[between(as.integer(group), 7, 9), age := fctr(3) ] - dt[between(as.integer(group), 10, 12), age := fctr(4) ] - dt[as.integer(group) >= 13, age := fctr(5) ] -} - cm_calc_R0_extended <- function( params ){ @@ -354,13 +343,11 @@ for(i in 1:nrow(run_options)){ model_seed = run_options[i, model_seed] )$dynamics - result <- reduce_ages( - sim[ + result <- sim[ compartment %in% c("cases","death_o","icu_p","nonicu_p","R") - ] - )[, + ][, .(value = sum(value)), - keyby = .(run, t, age, compartment) + keyby = .(run, t, group, compartment) ][value != 0][, run := i ] rm(params)