Merge pull request #35 from vimc/vimc-6748

Vimc 6748: vimpact flexibilities

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: vimpact
 Title: Vaccine Impact Calculation
-Version: 0.0.10
+Version: 0.1.0
 Authors@R:
     c(person(given = "Rich",
              family = "FitzJohn",

NEWS.md

@@ -1,3 +1,7 @@
+# vimpact 0.1.0
+
+* Increase flexibility - recognize new activity type, compare impact between different scenario types. 
+
 # vimpact 0.0.10
 
 * Include DTP3 in the disease vaccine delivery table to extract the coverage values from the DB

R/commonly_used_db_data.R

@@ -57,15 +57,20 @@ get_touchstone_id <- function(con, touchstone) {
 ##' 2. demographic uncertainty not only affects models, but also FVPs. If we are to conduct sensitivity analysis on impact_by_year_of_vaccination, we need to vary population input for adjusting FVPs.
 ##' @param full_description TRUE if including scenario_descriptions (coverage estimates will be duplicated for scenarios); and FALSE if only providing coverage estimates
 ##' @param demographic_source Demographic_source.code
+##' @param coverage_scenario_type Coverage scenario type. This is particularly useful for a coverage touchstone. Montagu coverage.name follows <disease>:<vaccine>,<activity_type>,<gavi_support_level>:<coverage_scenario_type>
+##' It is NULL by default. When it is not null, only pulls coverage.name that contains specified pattern.
 ##' @export
 extract_vaccination_history <- function(con, touchstone_cov = "201710gavi", touchstone_pop = NULL,
                                         year_min = 2000, year_max = 2100,
                                         vaccine_to_ignore = c("DTP3", "HepB_BD_home", "none"),
                                         disease_to_extract = NULL,
                                         countries_to_extract = NULL,
                                         gavi_support_levels = c("with", "bestminus"),
-                                        scenario_type = "default", external_population_estimates = NULL,
-                                        full_description = FALSE, demographic_source = NULL) {
+                                        scenario_type = "default",
+                                        external_population_estimates = NULL,
+                                        full_description = FALSE,
+                                        demographic_source = NULL,
+                                        coverage_scenario_type = NULL) {
 
   ## validate demography parameter
   ## when touchstone_pop is null, touchstone_cov is used for touchstone_pop
@@ -138,23 +143,25 @@ extract_vaccination_history <- function(con, touchstone_cov = "201710gavi", touc
                               touchstone_cov)
   if (nrow(cov_sets) == 0L) {
     # this is not a model run touchstone, need to extract coverage set directly
-    cov_sets2 <- DBI::dbGetQuery(con, "SELECT coverage_set.id AS coverage_set, vaccine, activity_type, gavi_support_level
+    cov_sets2 <- DBI::dbGetQuery(con, "SELECT name, coverage_set.id AS coverage_set, vaccine, activity_type, gavi_support_level
                                 FROM coverage_set
                                 WHERE touchstone = $1
                                 AND gavi_support_level != 'none'", touchstone_cov)
     cov_sets2 <- merge_by_common_cols(disease_vaccine_delivery, cov_sets2, all.y = TRUE)
-    cov_sets2$scenario_type <- "default"
+    if(!is.null(coverage_scenario_type)){
+      cov_sets2 <- cov_sets2[grepl(coverage_scenario_type, cov_sets2$name), ]
+    }
+    cov_sets2$scenario_type <- ifelse(is.null(coverage_scenario_type), "default", coverage_scenario_type)
     cov_sets2$scenario_description <- paste(cov_sets2$vaccine,
                                             cov_sets2$activity_type,
                                             cov_sets2$gavi_support_level,
                                             sep = "-")
-    cov_sets <- cov_sets2[names(cov_sets)]
+    cov_sets <- cov_sets2[c("scenario_type", "scenario_description", "disease", "coverage_set", "vaccine", "activity_type", "gavi_support_level")]
   }
 
   if (!is.null(disease_to_extract)) {
     cov_sets <- cov_sets[cov_sets$disease %in% disease_to_extract,]
   }
-
   country_ <- ifelse(is.null(countries_to_extract),
                      "",
                      sprintf("AND country IN %s", sql_in(countries_to_extract, text_item = TRUE)))
@@ -186,9 +193,10 @@ extract_vaccination_history <- function(con, touchstone_cov = "201710gavi", touc
   }
 
   ## spliting coverage data by age groups
-  cov1 <- cov[cov$activity_type == "routine",]
+  ll <- grepl("campaign", cov$activity_type)
+  cov1 <- cov[!ll,]
   cov1$age <- cov1$age_from
-  cov2 <- cov[cov$activity_type == "campaign",]
+  cov2 <- cov[ll,]
 
   d <- list(NULL)
   for (i in seq_along(cov2$activity_id)) {
@@ -229,7 +237,7 @@ extract_vaccination_history <- function(con, touchstone_cov = "201710gavi", touc
   cov2$country_nid <- country_id$nid[match(cov2$country, country_id$id)]
 
   if (!full_description) {
-    cov2$scenario_description <- "best-estimates"
+    cov2$scenario_description <- cov2$scenario_type
     cov2 <- unique(cov2)
   }
   cov2 <- cov2[order(cov2$scenario_description, cov2$delivery_id),]

R/impact_central.R

@@ -1,7 +1,7 @@
 get_raw_impact_details <- function(con, meta1, burden_outcome, is_under5 = FALSE, countries_to_extract = NULL){
   #verify parameters
   stopifnot(nrow(meta1) == 2L)
-  stopifnot(burden_outcome %in% c("deaths", "cases", "dalys"))
+  stopifnot(burden_outcome %in% c("deaths", "cases", "dalys", "yll"))
   stopifnot(is_under5 %in% c(TRUE, FALSE))
 
   #preparation
@@ -12,8 +12,10 @@ get_raw_impact_details <- function(con, meta1, burden_outcome, is_under5 = FALSE
     k <- 2
   } else if(burden_outcome == "dalys"){
     k <- 3
+  }else if(burden_outcome == "yll"){
+    k <- 4
   } else {
-    stop("Can only take burden outcome as one of deaths, cases, dalys")
+    stop("Can only take burden outcome as one of deaths, cases, dalys, yll")
   }
 
   # determine whether a recipe is for routine or campaign vaccine delivery
@@ -22,12 +24,15 @@ get_raw_impact_details <- function(con, meta1, burden_outcome, is_under5 = FALSE
   # if campaign, by calendar
   meta1$vaccine_delivery[meta1$vaccine_delivery == "no-vaccination"] <- ""
   v <- determine_vaccine_delivery(meta1)
-  i <- unique(grepl("routine", v))
-  j <- unique(grepl("campaign", v))
-  if(any(i) && any(j) && meta1$method[1] == "method2a"){
-    stop("method2a is vaccine delivery specific, cannot include both routine and campaign impact at the same time.")
+  if(length(v) > 0){
+    i <- unique(grepl("routine", v))
+    j <- unique(grepl("campaign", v))
+    if(any(i) && any(j) && meta1$method[1] == "method2a"){
+      stop("method2a is vaccine delivery specific, cannot include both routine and campaign impact at the same time.")
+    }
   }
 
+
   # constrain db extraction by country and age
 
   age_constrain <- ifelse(is_under5, "AND age < 5", "") #applies to all methods
@@ -136,7 +141,7 @@ impact_by_year_of_vaccination <- function(meta1, raw_impact, fvps, fvps_updates
   for (i in seq_along(v)){
     m <- unlist(strsplit(v[i], "-"))
     vaccine_delivery$vaccine[i] <- m[1]
-    vaccine_delivery$activity_type[i] <- m[2]
+    vaccine_delivery$activity_type[i] <- paste(m[-1], collapse = "-")
   }
   if (length(unique(vaccine_delivery$activity_type)) > 1L && method == "method2a"){
     stop("method2a can not accommodate both routine and campaign impact in the same time, as routine and campaign impact are calculated differently.")
@@ -212,10 +217,10 @@ determine_vaccine_delivery <- function(meta1){
 #' @return Impact ratio by country and burden outcome
 #' @export
 impact_by_year_of_vaccination_country_perspective <- function(
-  raw_impact, fvps, activity_type, vaccination_years) {
-  if (!(activity_type %in% c("routine", "campaign"))) {
+    raw_impact, fvps, activity_type, vaccination_years) {
+  if (!(activity_type %in% c("routine", "campaign", "routine-intensified"))) {
     stop(sprintf(
-      'Activity type must be "routine" or "campaign" got "%s".', activity_type))
+      'Activity type must be "routine", "campaign", "routine-intensified" got "%s".', activity_type))
   }
 
   ## Aggregate FVPs over years of vaccination
@@ -228,7 +233,7 @@ impact_by_year_of_vaccination_country_perspective <- function(
   ## Aggregate raw_impact grouped by country & burden_outcome where birth
   ## cohort is in range
   raw_impact$birth_cohort <- get_birth_cohort(raw_impact)
-  if (activity_type == "routine"){
+  if (grepl("routine", activity_type)){
     raw_impact <- raw_impact[raw_impact$birth_cohort %in%
                                (vaccination_years - min(fvps$age)), ]
   } else {
@@ -263,7 +268,7 @@ impact_by_year_of_vaccination_country_perspective <- function(
 #' @return Impact ratio by country, birth cohort and burden outcome
 #' @export
 impact_by_year_of_vaccination_cohort_perspective <- function(
-  raw_impact, fvps, vaccination_years) {
+    raw_impact, fvps, vaccination_years) {
 
   ## Aggregate FVPs by birth cohort and country
   fvps <- fvps[fvps$year %in% vaccination_years, ]
@@ -417,7 +422,7 @@ impact_by_birth_year <- function(baseline_burden, focal_burden) {
 #' and burden outcome
 #' @export
 impact_by_year_of_vaccination_activity_type <- function(
-  baseline_burden, focal_burden, fvps, vaccination_years) {
+    baseline_burden, focal_burden, fvps, vaccination_years) {
   assert_has_columns(
     baseline_burden,
     c("country", "burden_outcome", "year", "age", "value"))
@@ -451,7 +456,7 @@ impact_by_year_of_vaccination_activity_type <- function(
   }
 
   ## Routine
-  if (activity == "routine") {
+  if (grepl("routine", activity)) {
     raw_impact <- baseline_burden %>%
       impact_by_birth_year(focal_burden) %>%
       dplyr::rename(time = birth_cohort) %>%
@@ -515,7 +520,7 @@ impact_by_year_of_vaccination_activity_type <- function(
 #' and burden outcome
 #' @export
 impact_by_year_of_vaccination_birth_cohort <- function(
-  baseline_burden, focal_burden, fvps, vaccination_years) {
+    baseline_burden, focal_burden, fvps, vaccination_years) {
   assert_has_columns(
     baseline_burden,
     c("country", "burden_outcome", "year", "age", "value"))

R/impact_meta.R

@@ -61,7 +61,7 @@ get_meta_from_recipe <- function(default_recipe = TRUE, method = "method0", reci
   recipe$baseline <- gsub(";", ",", recipe$baseline)
   for (i in seq_along(recipe$burden_outcome)){
     if (recipe$burden_outcome[i] == "*"){
-      recipe$burden_outcome[i] <- "deaths;cases;dalys"
+      recipe$burden_outcome[i] <- "deaths;cases;dalys;yll"
     } else {
       recipe$burden_outcome[i] <- paste(recipe$burden_outcome[i], sep = ";")
     }
@@ -228,6 +228,11 @@ replace_burden_outcome <- function(burden_outcomes, a){
       m2[ii] <- m[k]
       ii <- ii+1
     }
+    k <- grepl("yll", v)
+    if(any(k)){
+      m2[ii] <- m[k]
+      ii <- ii+1
+    }
 
     t[i] <- paste(m2, collapse = ";")
   }

TODO.md

@@ -1,11 +1,3 @@
 ### TODOs
 
-* extract aggregated burden by recipe 
- - may be not necessary anymore, as this is output of method0 and method1 already, we can output it
- - but the output will only have counter-factual and best-estimate scenarios
- - if other scenarios are also needed, move burden calculation report here
-
-
-* standardise output to match past data sets - this is only for internal use, no need to export for non-science users
-
-* get_burden_central and get_burden_stochastic functions for easy-to-use
+Increase flexibility - allow addition activity types, comparison between different scenario types, external fvps for calculating impact ratios. 

inst/disease_vaccine_delivery.csv

@@ -1,24 +1,30 @@
-disease,vaccine,activity_type
-Cholera,Cholera,campaign
-HepB,HepB,routine
-HepB,HepB_BD,routine
-Hib,Hib3,routine
-HPV,HPV,campaign
-HPV,HPV,routine
-JE,JE,campaign
-JE,JE,routine
-Measles,MCV1,routine
-Measles,MCV2,routine
-Measles,Measles,campaign
-MenA,MenA,campaign
-MenA,MenA,routine
-PCV,PCV3,routine
-Rota,Rota,routine
-Rubella,RCV2,routine
-Rubella,Rubella,campaign
-Rubella,Rubella,routine
-Typhoid,Typhoid,campaign
-Typhoid,Typhoid,routine
-YF,YF,campaign
-YF,YF,routine
-DTP,DTP3,routine
+disease,vaccine,activity_type
+Cholera,Cholera,campaign
+DTP,DTP3,routine
+HepB,HepB,routine
+HepB,HepB_BD,routine
+HepB,HepB,routine-intensified
+Hib,Hib3,routine
+HPV,HPV,campaign
+HPV,HPV,routine
+HPV,HPV,routine-intensified
+JE,JE,campaign
+JE,JE,routine
+Measles,MCV1,routine
+Measles,MCV2,routine
+Measles,Measles,campaign
+Measles,Measles,routine-intensified
+MenA,MenA,campaign
+MenA,MenA,routine
+MenA,MenA,routine-intensified
+PCV,PCV3,routine
+Rota,Rota,routine
+Rubella,RCV2,routine
+Rubella,Rubella,campaign
+Rubella,Rubella,routine
+Rubella,Rubella,routine-intensified
+Typhoid,Typhoid,campaign
+Typhoid,Typhoid,routine
+YF,YF,campaign
+YF,YF,routine
+YF,YF,routine-intensified

tests/testthat/test-db-impact.R

@@ -175,7 +175,7 @@ test_that("impact calculation by year of vaccination country perspective", {
   expect_error(
     impact_by_year_of_vaccination_country_perspective(raw_impact, fvps, "test",
                                                       2000:2030),
-    'Activity type must be "routine" or "campaign" got "test".'
+    'Activity type must be "routine", "campaign", "routine-intensified" got "test".'
   )
 })