From ae15712b8d164257928f1a1ab9e4b99458859473 Mon Sep 17 00:00:00 2001
From: Wes Hinsley <w.hinsley@imperial.ac.uk>
Date: Thu, 19 Mar 2026 14:43:24 +0000
Subject: [PATCH 01/29] In progress

---
 DESCRIPTION          |  2 +-
 NEWS.md              |  9 ++++++--
 R/stochastic_files.R | 54 ++++++++++++++++++++++++++++++++++++++++----
 3 files changed, 57 insertions(+), 8 deletions(-)

diff --git a/DESCRIPTION b/DESCRIPTION
index 2e25561..1ab97ff 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: stoner
 Title: Support for Building VIMC Montagu Touchstones, using Dettl
-Version: 0.1.19
+Version: 0.1.20
 Authors@R:
     c(person("Wes", "Hinsley",role = c("aut", "cre", "cst", "dnc", "elg", "itr", "sng", "ard"),
            email = "w.hinsley@imperial.ac.uk"),
diff --git a/NEWS.md b/NEWS.md
index 0639fd8..86600b4 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,8 +1,13 @@
-# Stoner 0.1.19
+# Stoner 0.1.20
 
-* add stoner_stochastic_central to create central parquet from standardised stochastics.
+* Allow missing yll for processing older stochastics
+
+# Changes in 0.1.19
+
+* Add stoner_stochastic_central to create central parquet from standardised stochastics.
 
 # Changes in 0.1.18
+
 * Add stoner_stochastic_standardise for converting (drop)box stochastics to standard form
 * Add stoner_stochastic_graphs for quick plotting of stochastics.
 
diff --git a/R/stochastic_files.R b/R/stochastic_files.R
index 0f3f0db..922dda5 100644
--- a/R/stochastic_files.R
+++ b/R/stochastic_files.R
@@ -31,14 +31,19 @@
 ##' these to the simpler names. Processing Rubella stochastic files without
 ##' this set to TRUE will fail - so while we should always do this, keeping
 ##' the parameter makes it more clear in the code what we're doing and why.
+##' @param hepb2019_fix In 2019, HepB deaths and cases were subdivided into a
+##' number of different causes. This flag combines those into the single
+##' appropriate burden outcome.
 ##' @param missing_run_id_fix Some groups in the past have omitted run_id
 ##' from the files, but included them in the filenames. This fix inserts
 ##' that into the files if the index parameter indicates we have 200 runs to
 ##' process.
-
+##' @param allow_missing_yll yll was introduced in 2023? This flag allows
+##' it to be missing for processing older stochastics.
 stone_stochastic_standardise <- function(
     group, in_path, out_path, scenarios, files, index = 1,
-    rubella_fix = TRUE, missing_run_id_fix = TRUE) {
+    rubella_fix = TRUE, hepb2019_fix = TRUE, missing_run_id_fix = TRUE,
+    allow_missing_yll = TRUE) {
 
   dir.create(out_path, showWarnings = FALSE, recursive = TRUE)
   if ((length(files) == 1) && (grepl(":scenario", files))) {
@@ -56,19 +61,56 @@ stone_stochastic_standardise <- function(
       d <- read.csv(file.path(in_path, file))
       d$country_name <- NULL
 
-      # Fixes needed to standardise Rubella
+      # Various accumulated fixes for old/non-standard stochastics
+      # The flags are all true by default - Stoner will break untidily
+      # if the flags are turned off, but the problem occurs.
+
       if (rubella_fix) {
         names(d)[names(d) == "rubella_deaths_congenital"] <- "deaths"
         names(d)[names(d) == "rubella_cases_congenital"] <- "cases"
         d$rubella_infections <- NULL
       }
 
-      # Detect where run_id is missing, but in filenames
+      if (hepb2019_fix) {
+        if (!"cases" %in% names(d)) {
+          d$cases <- 0
+          cda_cases <- c("hepb_cases_acute_severe", "hepb_cases_dec_cirrh" ,
+                         "hepb_cases_hcc")
+          li_cases <- c("hepb_cases_acute_symp", "hepb_cases_fulminant",
+                        "hepb_cases_chronic",
+                        "hepb_chronic_symptomatic_in_acute_phase")
+          ic_cases <- c("hepb_cases_acute_severe", "hepb_cases_comp_cirrh",
+                        "hepb_cases_hcc_no_cirrh")
+          for (i in unique(c(cda_cases, li_cases, ic_cases))) {
+            if (i %in% names(d)) {
+              d$cases <- d$cases + d[[i]]
+              d[[i]] <- NULL
+            }
+          }
+        }
+
+        if (!"deaths" %in% names(d)) {
+          d$deaths <- 0
+          cda_deaths <- c("hepb_deaths_acute", "hepb_deaths_dec_cirrh",
+                          "hepb_deaths_hcc")
+          li_deaths <- c("hepb_deaths_acute", "hepb_deaths_total_cirrh",
+                         "hepb_deaths_hcc")
+          # IC just submitted deaths
+
+          for (i in unique(c(cda_deaths, li_deaths))) {
+            if (i %in% names(d)) {
+              d$deaths <- d$deaths + d[[i]]
+              d[[i]] <- NULL
+            }
+          }
+        }
+      }
 
       if (missing_run_id_fix) {
         if ((!"run_id" %in% names(d)) && (length(index) == 200)) d$run_id <- j
       }
 
+
       # Round to integer, as per guidance. (Not using as.integer, as that
       # has limits on how large numbers can be, so we are just truncating
       # digits here)
@@ -76,7 +118,9 @@ stone_stochastic_standardise <- function(
       d$dalys <- round(d$dalys)
       d$deaths <- round(d$deaths)
       d$cases <- round(d$cases)
-      d$yll <- round(d$yll)
+      if (("yll" %in% names(d)) || (!allow_missing_yll)) {
+        d$yll <- round(d$yll)
+      }
       d$cohort_size <- round(d$cohort_size)
 
       d

From 7c6d4ee282c85e19e4db3b430b3f5fa22c019a7b Mon Sep 17 00:00:00 2001
From: Wes Hinsley <w.hinsley@imperial.ac.uk>
Date: Tue, 24 Mar 2026 14:44:08 +0000
Subject: [PATCH 02/29] Fixes for early stochs

---
 R/stochastic_files.R | 70 +++++++++++++++++++++++++++++++++++++-------
 README.md            |  2 +-
 pkgdown/extra.css    |  4 +++
 3 files changed, 65 insertions(+), 11 deletions(-)

diff --git a/R/stochastic_files.R b/R/stochastic_files.R
index 922dda5..233851f 100644
--- a/R/stochastic_files.R
+++ b/R/stochastic_files.R
@@ -34,16 +34,25 @@
 ##' @param hepb2019_fix In 2019, HepB deaths and cases were subdivided into a
 ##' number of different causes. This flag combines those into the single
 ##' appropriate burden outcome.
+##' @param hib2019_fix In 2019, HepB deaths and cases were subdivided into a
+##' number of different causes. This flag combines those into the single
+##' appropriate burden outcome.
 ##' @param missing_run_id_fix Some groups in the past have omitted run_id
 ##' from the files, but included them in the filenames. This fix inserts
 ##' that into the files if the index parameter indicates we have 200 runs to
 ##' process.
 ##' @param allow_missing_yll yll was introduced in 2023? This flag allows
 ##' it to be missing for processing older stochastics.
+##' @param allow_missing_indexes In some early runs, different groups
+##' provided different numbers of files for each scenario, because some
+##' countries did not implement particular coverage campaigns. This
+##' flag needs to be TRUE for those groups, but the default is FALSE,
+##' since it's rare, and we generally want errors for missing files.
 stone_stochastic_standardise <- function(
     group, in_path, out_path, scenarios, files, index = 1,
-    rubella_fix = TRUE, hepb2019_fix = TRUE, missing_run_id_fix = TRUE,
-    allow_missing_yll = TRUE) {
+    rubella_fix = TRUE, hepb2019_fix = TRUE, hib2019_fix = TRUE,
+    missing_run_id_fix = TRUE, allow_missing_yll = TRUE,
+    allow_missing_dalys = TRUE, allow_missing_indexes = FALSE) {
 
   dir.create(out_path, showWarnings = FALSE, recursive = TRUE)
   if ((length(files) == 1) && (grepl(":scenario", files))) {
@@ -52,23 +61,35 @@ stone_stochastic_standardise <- function(
       files[j] <- gsub(":scenario", scenarios[j], files[j])
     }
   }
-
   for (i in seq_along(scenarios)) {
     message(scenarios[i])
     all_data <- as.data.frame(data.table::rbindlist(lapply(index, function(j) {
       cat("\r", j)
       file <- gsub(":index", j, files[i])
-      d <- read.csv(file.path(in_path, file))
+      filepath <- file.path(in_path, file)
+      if (!file.exists(filepath) && (allow_missing_indexes)) {
+        return(NULL)
+      }
+      d <- read.csv(filepath)
       d$country_name <- NULL
 
-      # Various accumulated fixes for old/non-standard stochastics
+      # Various accumulated fixes for e/non-standard stochastics
       # The flags are all true by default - Stoner will break untidily
       # if the flags are turned off, but the problem occurs.
 
       if (rubella_fix) {
-        names(d)[names(d) == "rubella_deaths_congenital"] <- "deaths"
-        names(d)[names(d) == "rubella_cases_congenital"] <- "cases"
-        d$rubella_infections <- NULL
+        if ("rubella_deaths_congenital" %in% names(d)) {
+          message("Converting rubella_deaths_congenital to deaths")
+          names(d)[names(d) == "rubella_deaths_congenital"] <- "deaths"
+        }
+        if ("rubella_cases_congenital" %in% names(d)) {
+          message("Converting rubella_cases_congenital to cases")
+          names(d)[names(d) == "rubella_cases_congenital"] <- "cases"
+        }
+        if ("rubella_infections" %in% names(d)) {
+          message("Ignoring rubella_infections")
+          d$rubella_infections <- NULL
+        }
       }
 
       if (hepb2019_fix) {
@@ -83,6 +104,7 @@ stone_stochastic_standardise <- function(
                         "hepb_cases_hcc_no_cirrh")
           for (i in unique(c(cda_cases, li_cases, ic_cases))) {
             if (i %in% names(d)) {
+              message(sprintf("Including %s in cases", i))
               d$cases <- d$cases + d[[i]]
               d[[i]] <- NULL
             }
@@ -95,10 +117,12 @@ stone_stochastic_standardise <- function(
                           "hepb_deaths_hcc")
           li_deaths <- c("hepb_deaths_acute", "hepb_deaths_total_cirrh",
                          "hepb_deaths_hcc")
-          # IC just submitted deaths
+          ic_deaths <- c("hepb_deaths_acute", "hepb_deaths_comp_cirrh",
+                         "hepb_deaths_dec_cirrh", "hepb_deaths_hcc")
 
           for (i in unique(c(cda_deaths, li_deaths))) {
             if (i %in% names(d)) {
+              message(sprintf("Including %s in deaths", i))
               d$deaths <- d$deaths + d[[i]]
               d[[i]] <- NULL
             }
@@ -106,6 +130,25 @@ stone_stochastic_standardise <- function(
         }
       }
 
+      if (hib2019_fix) {
+        if (("cases_pneumo" %in% names(d)) &&
+            ("cases_men" %in% names(d)) &&
+            (!"cases" %in% names(d))) {
+          message("cases = cases_men + cases_pneumo")
+          d$cases <- d$cases_pneumo + d$cases_men
+          d$cases_pneumo <- NULL
+          d$cases_men <- NULL
+        }
+        if (("deaths_pneumo" %in% names(d)) &&
+            ("deaths_men" %in% names(d)) &&
+            (!"deaths" %in% names(d))) {
+          message("deaths = deaths_men + deaths_pneumo")
+          d$deaths <- d$deaths_pneumo + d$deaths_men
+          d$deaths_pneumo <- NULL
+          d$deaths_men <- NULL
+        }
+      }
+
       if (missing_run_id_fix) {
         if ((!"run_id" %in% names(d)) && (length(index) == 200)) d$run_id <- j
       }
@@ -115,11 +158,18 @@ stone_stochastic_standardise <- function(
       # has limits on how large numbers can be, so we are just truncating
       # digits here)
 
-      d$dalys <- round(d$dalys)
+      if (("dalys" %in% names(d)) || (!allow_missing_dalys)) {
+        d$dalys <- round(d$dalys)
+      } else {
+        message("Dalys missing. (Ignored)")
+      }
+
       d$deaths <- round(d$deaths)
       d$cases <- round(d$cases)
       if (("yll" %in% names(d)) || (!allow_missing_yll)) {
         d$yll <- round(d$yll)
+      } else {
+        message("yll missing. (Ignored)")
       }
       d$cohort_size <- round(d$cohort_size)
 
diff --git a/README.md b/README.md
index c115f9f..163a6e3 100644
--- a/README.md
+++ b/README.md
@@ -4,7 +4,7 @@
 [![R build status](https://github.com/vimc/stoner/workflows/R-CMD-check/badge.svg)](https://github.com/vimc/stoner/actions)
 <!-- badges: end -->
 
-<pre>
+<pre style="line-height: 1;">
      .-'''-. ,---------.    ,-----.    ,---.   .--.    .-''-.  .-------.
     / _     \\          \ .'  .-,  '.  |    \  |  |  .'_ _   \ |  _ _   \
    (`' )/`--' `--.  ,---'/ ,-.|  \ _ \ |  ,  \ |  | / ( ` )   '| ( ' )  |
diff --git a/pkgdown/extra.css b/pkgdown/extra.css
index a5229e3..4a0ae85 100644
--- a/pkgdown/extra.css
+++ b/pkgdown/extra.css
@@ -1,3 +1,7 @@
 .row > main {
   hyphens: none;
 }
+
+pre {
+  line-height: 1;
+}

From 8de250d554288b2924c22e194303c65912f29fa6 Mon Sep 17 00:00:00 2001
From: Wes Hinsley <w.hinsley@imperial.ac.uk>
Date: Thu, 26 Mar 2026 11:23:01 +0000
Subject: [PATCH 03/29] Fix zero burdens with LiST PCV/Hib

---
 R/stochastic_files.R | 7 +++----
 1 file changed, 3 insertions(+), 4 deletions(-)

diff --git a/R/stochastic_files.R b/R/stochastic_files.R
index 233851f..0ed805c 100644
--- a/R/stochastic_files.R
+++ b/R/stochastic_files.R
@@ -132,16 +132,15 @@ stone_stochastic_standardise <- function(
 
       if (hib2019_fix) {
         if (("cases_pneumo" %in% names(d)) &&
-            ("cases_men" %in% names(d)) &&
-            (!"cases" %in% names(d))) {
+            ("cases_men" %in% names(d))) {
+
           message("cases = cases_men + cases_pneumo")
           d$cases <- d$cases_pneumo + d$cases_men
           d$cases_pneumo <- NULL
           d$cases_men <- NULL
         }
         if (("deaths_pneumo" %in% names(d)) &&
-            ("deaths_men" %in% names(d)) &&
-            (!"deaths" %in% names(d))) {
+            ("deaths_men" %in% names(d))) {
           message("deaths = deaths_men + deaths_pneumo")
           d$deaths <- d$deaths_pneumo + d$deaths_men
           d$deaths_pneumo <- NULL

From 58d589a760d130b7df198f9750d22f878b49b997 Mon Sep 17 00:00:00 2001
From: Wes Hinsley <w.hinsley@imperial.ac.uk>
Date: Thu, 26 Mar 2026 11:23:44 +0000
Subject: [PATCH 04/29] Stochastic scripts 2017-2023

---
 scripts/process_stoch_201710gavi.R | 298 ++++++++++++++++++++
 scripts/process_stoch_201910gavi.R | 419 +++++++++++++++++++++++++++++
 scripts/process_stoch_202110gavi.R | 386 ++++++++++++++++++++++++++
 scripts/process_stoch_202310gavi.R |   2 +-
 4 files changed, 1104 insertions(+), 1 deletion(-)
 create mode 100644 scripts/process_stoch_201710gavi.R
 create mode 100644 scripts/process_stoch_201910gavi.R
 create mode 100644 scripts/process_stoch_202110gavi.R

diff --git a/scripts/process_stoch_201710gavi.R b/scripts/process_stoch_201710gavi.R
new file mode 100644
index 0000000..bcb89f4
--- /dev/null
+++ b/scripts/process_stoch_201710gavi.R
@@ -0,0 +1,298 @@
+base_in_path <- "//wpia-hn2.hpc.dide.ic.ac.uk/vimc_stochastics_dropbox/latest/201710gavi"
+base_out_path <- "//wpia-hn2.hpc.dide.ic.ac.uk/vimc_stochastics/201710gavi"
+
+# In case useful for looking up scenarios
+
+vault <- vaultr::vault_client(login = "github")
+password <- vault$read("/secret/vimc/database/production/users/readonly")$password
+con <- DBI::dbConnect(RPostgres::Postgres(),
+                      dbname = "montagu",
+                      host = "montagu.vaccineimpact.org",
+                      port = 5432, password = password,
+                      user = "readonly")
+
+fetch_scenarios <- function(disease, touchstone='201710gavi-5') {
+  sort(DBI::dbGetQuery(con, "
+    SELECT scenario_description FROM
+           scenario JOIN scenario_description
+           ON scenario.scenario_description = scenario_description.id
+           WHERE disease = $1 AND touchstone = $2", list(disease, touchstone))$scenario_description)
+}
+
+# Let's unleash the cow
+
+setwd("Q:/testcow")
+writeLines("vimc/stoner@VIMC-9230", "pkgdepends.txt")
+hipercow::hipercow_init(driver = "dide-windows")
+hipercow::hipercow_provision()
+# Network/memory might be too much for more than a job per node.
+hres <- hipercow::hipercow_resources(cores = 1L, exclusive = TRUE)
+
+
+###############
+# HepB
+
+scenarios <- c("hepb-no-vaccination",
+               "hepb-hepb-routine-with",
+               "hepb-bd-routine-with",
+               "hepb-bd-routine-with-hepb-routine-with",
+               "hepb-bd-routine-best-hepb-routine-with")
+
+stoner::stone_stochastic_standardise(
+  group = "CDA-Razavi",
+  in_path = file.path(base_in_path, "CDA-Razavi"),
+  out_path = file.path(base_out_path, "HepB_CDA-Razavi"),
+  scenarios = scenarios,
+  files = "Devin Razavi-Shearer - stochastic_burden_template_HepB-CDA-Razavi_:scenario_:index.csv.xz",
+  index = 1:98,
+  allow_missing_indexes = TRUE)
+
+hipercow::task_create_expr(resources = hres, expr =
+  stoner::stone_stochastic_standardise(
+    group = "Li",
+    in_path = file.path(base_in_path, "Li"),
+    out_path = file.path(base_out_path, "HepB_Li"),
+    scenarios = scenarios,
+    files = ":scenario:index.csv.xz",
+    index = 1:200))
+
+hipercow::task_create_expr(resources = hres, expr =
+  stoner::stone_stochastic_standardise(
+    group = "IC-Hallett",
+    in_path = file.path(base_in_path, "IC-Hallett"),
+    out_path = file.path(base_out_path, "HepB_IC-Hallett"),
+    scenarios = scenarios,
+    files = "stochastic_burden_est_:scenario_:index.csv.xz",
+    index = 1:200))
+
+
+###############
+# HIB
+
+scenarios <- c("hib-no-vaccination", "hib-routine-gavi")
+
+stoner::stone_stochastic_standardise(
+  group = "JHU-Tam",
+  in_path = file.path(base_in_path, "JHU-Tam-Hib"),
+  out_path = file.path(base_out_path, "Hib_JHU-Tam"),
+  scenarios = scenarios,
+  files = c("No_Hib_all.csv.xz", "With_Hib_all.csv.xz"))
+
+stoner::stone_stochastic_standardise(
+  group = "LSHTM-Clark",
+  in_path = file.path(base_in_path, "LSHTM-Clark-Hib"),
+  out_path = file.path(base_out_path, "Hib_LSHTM-Clark"),
+  scenarios = scenarios,
+  files = c("Kaja Abbas - stochastic-burden-estimate.201710gavi-5.Hib_LSHTM-CLark_standard-Hib-no-vaccination.csv.xz",
+            "Kaja Abbas - stochastic-burden-estimate.201710gavi-5.Hib_LSHTM-CLark_standard-Hib-routine-gavi.csv.xz"))
+
+
+###############
+# HPV
+
+scenarios <- c("hpv-no-vaccination", "hpv-routine-gavi",
+               "hpv-campaign-gavi")
+
+hipercow::task_create_expr(resources = hres, expr =
+stoner::stone_stochastic_standardise(
+  group = "Harvard-Sweet",
+  in_path = file.path(base_in_path, "Harvard-Sweet"),
+  out_path = file.path(base_out_path, "HPV_Harvard-Sweet"),
+  scenarios = scenarios,
+  files = c("1. central_burden_template_HPV-Harvard-Sweet No Vaccine_PSA-:index_08.22.19.csv.xz",
+            "2. central_burden_template_HPV-Harvard-Sweet Routine_PSA-:index_09.09.19.csv.xz",
+            "3. central_burden_template_HPV-Harvard-Sweet Campaign_PSA-:index_08.22.19.csv.xz"),
+  index = 1:200))
+
+hipercow::task_create_expr(resources = hres, expr =
+stoner::stone_stochastic_standardise(
+  group = "LSHTM-Jit",
+  in_path = file.path(base_in_path, "LSHTM-Jit-HPV"),
+  out_path = file.path(base_out_path, "HPV_LSHTM-Jit"),
+  scenarios =  scenarios,
+  files = "Kaja Abbas - stochastic-burden-estimate.201710gavi-5.HPV_LSHTM-Jit_standard-:scenario.csv.xz"))
+
+###############
+# JE
+
+scenarios <- c("je-routine-no-vaccination", "je-campaign-gavi",
+               "je-routine-gavi")
+
+hipercow::task_create_expr(resources = hres, expr =
+stoner::stone_stochastic_standardise(
+  group = "OUCRU-Clapham",
+  in_path = file.path(base_in_path, "OUCRU-Clapham"),
+  out_path = file.path(base_out_path, "JE_OUCRU-Clapham"),
+  scenarios =  scenarios,
+  files = "Duy Nguyen - stochastic-burden-template.201710gavi-6.JE_OUCRU-Clapham_standard_:scenario_:index.csv.xz",
+  index = 1:200))
+
+stub <- "Sean Moore - stochastic_burden_est_JE_UND-Moore"
+
+hipercow::task_create_expr(resources = hres, expr =
+stoner::stone_stochastic_standardise(
+  group = "UND-Moore",
+  in_path = file.path(base_in_path, "UND-Moore"),
+  out_path = file.path(base_out_path, "JE_UND-Moore"),
+  scenarios =  scenarios,
+  files = c("Sean Moore - stochastic_burden_est_JE_UND-Moore_je-no-vaccination.csv.xz",
+            "Sean Moore - stochastic_burden_est_JE_UND-Moore_je-campaign-gavi.csv.xz",
+            "Sean Moore - stochastic_burden_est_JE_UND-Moore_je-routine-gavi.csv.xz")))
+
+
+
+###############
+# Measles
+scenarios <- c("measles-no-vaccination", "measles-mcv1-gavi",
+               "measles-mcv2-gavi", "measles-campaign-gavi")
+
+hipercow::task_create_expr(resources = hres, expr =
+stoner::stone_stochastic_standardise(
+  group = "PSU-Ferrari",
+  in_path = file.path(base_in_path, "PSU-Ferrari"),
+  out_path = file.path(base_out_path, "Measles_PSU-Ferrari"),
+  scenarios = scenarios,
+  files = c("Matthew Ferrari - stochastic_burden_novax_Measles-PSU-Ferrari-5.csv.xz",
+            "Matthew Ferrari - stochastic_burden_mcv1_Measles-PSU-Ferrari-5.csv.xz",
+            "Matthew Ferrari - stochastic_burden_mcv2_Measles-PSU-Ferrari-5.csv.xz",
+            "Matthew Ferrari - stochastic_burden_sia_Measles-PSU-Ferrari-5.csv.xz")))
+
+hipercow::task_create_expr(resources = hres, expr =
+stoner::stone_stochastic_standardise(
+  group = "LSHTM-Jit",
+  in_path = file.path(base_in_path, "LSHTM-Jit-Measles"),
+  out_path = file.path(base_out_path, "Measles_LSHTM-Jit"),
+  scenarios = scenarios,
+  files = c("stochastic_burden_template_Measles-LSHTM-Jit-no_vaxx.csv.xz",
+            "Petra Klepac - stochastic_burden_template_Measles-LSHTM-Jit-mcv1-only.csv.xz",
+            "stochastic_burden_template_Measles-LSHTM-Jit-mcv1_mcv2.csv.xz",
+            "stochastic_burden_template_Measles-LSHTM-Jit-mcv1-mcv2-campaigns.csv.xz")))
+
+###############
+# MenA
+
+scenarios = c("mena-no-vaccination", "mena-campaign-gavi",
+              "mena-routine-gavi")
+
+stoner::stone_stochastic_standardise(
+  group = "Cambridge-Trotter",
+  in_path = file.path(base_in_path, "Cambridge-Trotter"),
+  out_path = file.path(base_out_path, "MenA_Cambridge-Trotter"),
+  scenarios = scenarios,
+  files = c("Cambridge-Trotter-mena-no-vacc.csv.xz",
+            "Cambridge-Trotter-mena-campaign.csv.xz",
+            "Cambridge-Trotter-mena-routine.csv.xz"))
+
+stoner::stone_stochastic_standardise(
+  group = "KPW-Jackson",
+  in_path = file.path(base_in_path, "KPW-Jackson"),
+  out_path = file.path(base_out_path, "MenA_KPW-Jackson"),
+  scenarios = scenarios,
+  files = c("Mike Jackson - stochastic_burden_est_MenA-KPW-Jackson_mena-no-vaccination_:index.csv.xz",
+            "Mike Jackson - stochastic_burden_est_MenA-KPW-Jackson_mena-campaign-gavi_:index.csv.xz",
+            "Mike Jackson - stochastic_burden_est_MenA-KPW-Jackson_mena-routine-gavi_:index.csv.xz"),
+  index = 1:2,
+  allow_missing_indexes = TRUE)
+
+###############
+# PCV
+
+scenarios <- c("pcv-no-vaccination", "pcv-routine-gavi")
+
+stoner::stone_stochastic_standardise(
+  group = "JHU-Tam",
+  in_path = file.path(base_in_path, "JHU-Tam-PCV"),
+  out_path = file.path(base_out_path, "PCV_JHU-Tam"),
+  scenarios = scenarios,
+  files = c("No_PCV_all.csv.xz", "With_PCV_all.csv.xz"))
+
+stoner::stone_stochastic_standardise(
+  group = "LSHTM-Clark",
+  in_path = file.path(base_in_path, "LSHTM-Clark-PCV"),
+  out_path = file.path(base_out_path, "PCV_LSHTM-Clark"),
+  scenarios = scenarios,
+  files = c("Kaja Abbas - stochastic-burden-estimate.201710gavi-5.PCV_LSHTM-CLark_standard-PCV-no-vaccination.csv.xz",
+            "Kaja Abbas - stochastic-burden-estimate.201710gavi-5.PCV_LSHTM-CLark_standard-PCV-routine-gavi.csv.xz"))
+
+###############
+# Rota
+
+scenarios <- c("rota-no-vaccination", "rota-routine-gavi")
+
+stoner::stone_stochastic_standardise(
+  group = "JHU-Tam",
+  in_path = file.path(base_in_path, "JHU-Tam-Rota"),
+  out_path = file.path(base_out_path, "Rota_JHU-Tam"),
+  scenarios = scenarios,
+  files = c("No_Rota_all.csv.xz", "With_Rota_all.csv.xz"))
+
+stoner::stone_stochastic_standardise(
+  group = "LSHTM-Clark",
+  in_path = file.path(base_in_path, "LSHTM-Clark-Rota"),
+  out_path = file.path(base_out_path, "Rota_LSHTM-Clark"),
+  scenarios = scenarios,
+  files = c("Kaja Abbas - stochastic-burden-estimate.201710gavi-5.Rota_LSHTM-CLark_standard-Rota-no-vaccination.csv.xz",
+            "Kaja Abbas - stochastic-burden-estimate.201710gavi-5.Rota_LSHTM-CLark_standard-Rota-routine-gavi.csv.xz"))
+
+###############
+# Rubella
+
+scenarios <- c(
+  "rubella-routine-no-vaccination", "rubella-campaign-gavi",
+  "rubella-rcv1-gavi", "rubella-rcv2-gavi")
+
+#hipercow::task_create_expr(resources = hres, expr =
+stoner::stone_stochastic_standardise(
+  group = "JHU-Lessler",
+  in_path = file.path(base_in_path, "JHU-Lessler"),
+  out_path = file.path(base_out_path, "Rubella_JHU-Lessler"),
+  scenarios = scenarios,
+  files = c("Amy Winter - stochastic_burden_est-rubella-no-vaccination_:index.csv.xz",
+            "Amy Winter - stochastic_burden_est-rubella-campaign-gavi_:index.csv.xz",
+            "Amy Winter - stochastic_burden_est-rubella-routine-gavi_:index.csv.xz",
+            "Amy Winter - stochastic_burden_est-rubella-rcv2-gavi_:index.csv.xz"),
+  index = 1:10)
+
+scenarios <- c(
+  "rubella-routine-no-vaccination", "rubella-campaign-gavi",
+  "rubella-rcv2-gavi")
+
+hipercow::task_create_expr(resources = hres, expr =
+stoner::stone_stochastic_standardise(
+  group = "PHE-Vynnycky",
+  in_path = file.path(base_in_path, "PHE-Vynnycky"),
+  out_path = file.path(base_out_path, "Rubella_PHE-Vynnycky"),
+  scenarios = scenarios,
+  files = c("stochastic_burden_estimate_Rubella-PHE-Vynnycky_rubella-routine-no-vaccination_:index.csv.xz",
+            "stochastic_burden_estimate_Rubella-PHE-Vynnycky_rubella-campaign-gavi_:index.csv.xz",
+            "stochastic_burden_estimate_Rubella-PHE-Vynnycky_rubella-routine-gavi_:index.csv.xz"),
+  index = 1:200))
+
+
+###############
+# YF
+
+scenarios <- c("yf-no-vaccination", "yf-preventive-gavi",
+               "yf-routine-gavi")
+
+hipercow::task_create_expr(resources = hres, expr =
+stoner::stone_stochastic_standardise(
+  group = "IC-Garske",
+  in_path = file.path(base_in_path, "IC-Garske"),
+  out_path = file.path(base_out_path, "YF_IC-Garske"),
+  scenarios = scenarios,
+  files = c("Katy Gaythorpe - NEW_NEW_stochastic_burden_estimate_YF_no-vaccination_:index.csv.xz",
+            "Katy Gaythorpe - NEW_NEW_stochastic_burden_estimate_YF_preventive-gavi_:index.csv.xz",
+            "Katy Gaythorpe - NEW_NEW_stochastic_burden_estimate_YF_routine-gavi_:index.csv.xz"),
+  index = 1:200))
+
+hipercow::task_create_expr(resources = hres, expr =
+stoner::stone_stochastic_standardise(
+  group = "UND-Perkins",
+  in_path = file.path(base_in_path, "UND-Perkins"),
+  out_path = file.path(base_out_path, "YF_UND-Perkins"),
+  scenarios = scenarios,
+  files = c("John Huber - stochastic_burden_est_YF_UND-Perkins_yf-no-vaccination-gavi_:index.csv.xz",
+            "John Huber - stochastic_burden_est_YF_UND-Perkins_yf-preventive-gavi_:index.csv.xz",
+            "John Huber - stochastic_burden_est_YF_UND-Perkins_yf-routine-gavi_:index.csv.xz"),
+  index = 1:200))
diff --git a/scripts/process_stoch_201910gavi.R b/scripts/process_stoch_201910gavi.R
new file mode 100644
index 0000000..1ed1318
--- /dev/null
+++ b/scripts/process_stoch_201910gavi.R
@@ -0,0 +1,419 @@
+base_in_path <- "//wpia-hn2.hpc.dide.ic.ac.uk/vimc_stochastics_dropbox/latest/201910gavi"
+base_out_path <- "//wpia-hn2.hpc.dide.ic.ac.uk/vimc_stochastics/201910gavi"
+
+# In case useful for looking up scenarios
+
+vault <- vaultr::vault_client(login = "github")
+password <- vault$read("/secret/vimc/database/production/users/readonly")$password
+con <- DBI::dbConnect(RPostgres::Postgres(),
+                      dbname = "montagu",
+                      host = "montagu.vaccineimpact.org",
+                      port = 5432, password = password,
+                      user = "readonly")
+
+fetch_scenarios <- function(disease) {
+  sort(DBI::dbGetQuery(con, "
+    SELECT scenario_description FROM
+           scenario JOIN scenario_description
+           ON scenario.scenario_description = scenario_description.id
+           WHERE disease = $1 AND touchstone='201910gavi-5'", disease)$scenario_description)
+}
+
+# Let's unleash the cow
+
+setwd("Q:/testcow")
+writeLines("vimc/stoner@VIMC-9230", "pkgdepends.txt")
+hipercow::hipercow_init(driver = "dide-windows")
+hipercow::hipercow_provision()
+# Network/memory might be too much for more than a job per node.
+hres <- hipercow::hipercow_resources(cores = 1L, exclusive = TRUE)
+
+###############
+# Cholera
+
+scenarios = c("cholera-no-vaccination", "cholera-campaign-default")
+
+stoner::stone_stochastic_standardise(
+  group = "IVI-Kim",
+  in_path = file.path(base_in_path, "IVI-Kim-Cholera"),
+  out_path = file.path(base_out_path, "Cholera_IVI-Kim"),
+  scenarios = scenarios,
+  files = c("Jong-Hoon Kim - stoch_output_Cholera_novacc_20210902.csv.xz",
+            "Jong-Hoon Kim - stoch_output_Cholera_campaign_20210902.csv.xz")
+  )
+
+stoner::stone_stochastic_standardise(
+  group = "JHU-Lee",
+  in_path = file.path(base_in_path, "JHU-Lee"),
+  out_path = file.path(base_out_path, "Cholera_JHU-Lee"),
+  scenarios = scenarios,
+  files = c("Kaiyue Zou - stochastic-burden-template.201910gavi-5.Cholera_no-vaccination.csv.xz",
+            "Kaiyue Zou - stochastic-burden-template.201910gavi-5.Cholera_campaign-default.csv.xz")
+)
+
+###############
+# HepB
+
+scenarios <- c("hepb-no-vaccination",
+               "hepb-bd-default-hepb-routine-default",
+               "hepb-bd-routine-bestcase-hepb-routine-bestcase",
+               "hepb-bd-routine-bestcase",
+               "hepb-bd-routine-default",
+               "hepb-hepb-routine-bestcase",
+               "hepb-hepb-routine-default",
+               "hepb-stop")
+
+stub <- "Ivane Gamkrelidze - stochastic-burden-template.201910gavi-4.HepB_CDA-Razavi"
+
+stoner::stone_stochastic_standardise(
+  group = "CDA-Razavi",
+  in_path = file.path(base_in_path, "CDA-Razavi"),
+  out_path = file.path(base_out_path, "HepB_CDA-Razavi"),
+  scenarios = scenarios,
+  files =
+    c(sprintf("%s_all_hepb-no-vaccination.csv.xz", stub),
+      sprintf("%s_all_hepb-bd-default-hepb-routine-default.csv.xz", stub),
+      sprintf("%s_all_hepb-bd-routine-bestcase-hepb-routine-bestcase.csv.xz", stub),
+      sprintf("%s_bd_hepb-bd-routine-bestcase.csv.xz", stub),
+      sprintf("%s_bd_hepb-bd-routine-default.csv.xz", stub),
+      sprintf("%s_non_bd_hepb-hepb-routine-bestcase.csv.xz", stub),
+      sprintf("%s_non_bd_hepb-hepb-routine-default.csv.xz", stub),
+      sprintf("%s_all_hepb-stop.csv.xz", stub))
+)
+
+hipercow::task_create_expr(resources = hres, expr =
+  stoner::stone_stochastic_standardise(
+    group = "Li",
+    in_path = file.path(base_in_path, "Li"),
+    out_path = file.path(base_out_path, "HepB_Li"),
+    scenarios = scenarios,
+    files = ":scenario:index.csv.xz",
+    index = 1:200))
+
+hipercow::task_create_expr(resources = hres, expr =
+  stoner::stone_stochastic_standardise(
+    group = "IC-Hallett",
+    in_path = file.path(base_in_path, "IC-Hallett"),
+    out_path = file.path(base_out_path, "HepB_IC-Hallett"),
+    scenarios = scenarios,
+    files = "stochastic_burden_est_HepB-IC-Hallett_:scenario_:index.csv.xz",
+    index = 1:200))
+
+
+###############
+# HIB
+
+scenarios <- c("hib-no-vaccination", "hib-routine-default",
+               "hib-routine-bestcase")
+
+stoner::stone_stochastic_standardise(
+  group = "JHU-Tam",
+  in_path = file.path(base_in_path, "JHU-Tam-Hib"),
+  out_path = file.path(base_out_path, "Hib_JHU-Tam"),
+  scenarios = scenarios,
+  files = c("novac:index.csv.xz", "default:index.csv.xz", "best:index.csv.xz"),
+  index = 1:14)
+
+hipercow::task_create_expr(resources = hres, expr =
+  stoner::stone_stochastic_standardise(
+  group = "LSHTM-Clark",
+  in_path = file.path(base_in_path, "LSHTM-Clark_Hib"),
+  out_path = file.path(base_out_path, "Hib_LSHTM-Clark"),
+  scenarios = scenarios,
+  files = c("VIMC_Hib_PSA_NoVax.csv.xz", "VIMC_Hib_PSA_Default.csv.xz",
+            "VIMC_Hib_PSA_Best.csv.xz")))
+
+###############
+# HPV
+
+scenarios <- c("hpv-no-vaccination", "hpv-routine-bestcase",
+               "hpv-campaign-bestcase", "hpv-routine-default",
+               "hpv-campaign-default")
+
+stub <- "stochastic-burden-est.201910gavi-5.HPV_Harvard-Sweet"
+
+stoner::stone_stochastic_standardise(
+  group = "Harvard-Sweet",
+  in_path = file.path(base_in_path, "Harvard-Sweet"),
+  out_path = file.path(base_out_path, "HPV_Harvard-Sweet"),
+  scenarios = scenarios,
+  files = c(sprintf("%s_novacc_run_:index.csv.xz", stub),
+            sprintf("%s_routine-bestcase_run_:index.csv.xz", stub),
+            sprintf("%s_campaign-bestcase_run_:index.csv.xz", stub),
+            sprintf("%s_routine-default_run_:index.csv.xz", stub),
+            sprintf("%s_campaign-default_run_:index.csv.xz", stub)),
+  index = 1:200)
+
+hipercow::task_create_expr(resources = hres, expr =
+stoner::stone_stochastic_standardise(
+  group = "LSHTM-Jit",
+  in_path = file.path(base_in_path, "LSHTM-Jit_HPV"),
+  out_path = file.path(base_out_path, "HPV_LSHTM-Jit"),
+  scenarios =  scenarios,
+  files =
+    c("stochastic-burden-novaccination_201910gavi-4_hpv-no-vaccination.csv.xz",
+      "stochastic-burden-vaccination_201910gavi-4_hpv-campaign-bestcase.csv.xz",
+      "stochastic-burden-vaccination_201910gavi-4_hpv-campaign-default.csv.xz",
+      "stochastic-burden-vaccination_201910gavi-4_hpv-routine-bestcase.csv.xz",
+      "stochastic-burden-vaccination_201910gavi-4_hpv-routine-default.csv.xz")))
+
+###############
+# JE
+
+scenarios <- c("je-routine-no-vaccination", "je-campaign-default",
+               "je-campaign-bestcase", "je-routine-default",
+               "je-routine-bestcase")
+
+hipercow::task_create_expr(resources = hres, expr =
+stoner::stone_stochastic_standardise(
+  group = "OUCRU-Clapham",
+  in_path = file.path(base_in_path, "OUCRU-Clapham"),
+  out_path = file.path(base_out_path, "JE_OUCRU-Clapham"),
+  scenarios =  scenarios,
+  files = c(
+    "Template_Stochastic_Naive4_correcting_:index.csv.xz",
+    "Template_Stochastic_Campaign_Default4_correcting_:index.csv.xz",
+    "Template_Stochastic_Campaign_Best4_correcting_:index.csv.xz",
+    "Template_Stochastic_Routine_Default4_correcting_:index.csv.xz",
+    "Template_Stochastic_Routine_Best4_correcting_:index.csv.xz"),
+  index = 1:200))
+
+stub <- "Sean Moore - stochastic_burden_est_JE_UND-Moore"
+
+hipercow::task_create_expr(resources = hres, expr =
+stoner::stone_stochastic_standardise(
+  group = "UND-Moore",
+  in_path = file.path(base_in_path, "UND-Moore"),
+  out_path = file.path(base_out_path, "JE_UND-Moore"),
+  scenarios =  scenarios,
+  files = c(
+    sprintf("%s_je-campaign-bestcase.csv.xz", stub),
+    sprintf("%s_je-campaign-default.csv.xz", stub),
+    sprintf("%s_je-no-vaccination.csv.xz", stub),
+    sprintf("%s_je-routine-bestcase.csv.xz", stub),
+    sprintf("%s_je-routine-default.csv.xz", stub))))
+
+
+
+###############
+# Measles
+scenarios <- c("measles-no-vaccination", "measles-mcv1-default",
+               "measles-mcv2-default", "measles-mcv1-bestcase",
+               "measles-mcv2-bestcase", "measles-campaign-default",
+               "measles-campaign-only-default", "measles-campaign-bestcase",
+               "measles-campaign-only-bestcase", "measles-stop")
+
+hipercow::task_create_expr(resources = hres, expr =
+stoner::stone_stochastic_standardise(
+  group = "PSU-Ferrari",
+  in_path = file.path(base_in_path, "PSU-Ferrari"),
+  out_path = file.path(base_out_path, "Measles_PSU-Ferrari"),
+  scenarios = scenarios,
+  files = c(
+    "Heather Santos - novax_stochastic:index_burden_Measles-PSU-Ferrari.csv.xz",
+    "Heather Santos - default_mcv1_stochastic:index_burden_Measles-PSU-Ferrari.csv.xz",
+    "Heather Santos - default_mcv2_stochastic:index_burden_Measles-PSU-Ferrari.csv.xz",
+    "Heather Santos - bestcase_mcv1_stochastic:index_burden_Measles-PSU-Ferrari.csv.xz",
+    "Heather Santos - bestcase_mcv2_stochastic:index_burden_Measles-PSU-Ferrari.csv.xz",
+    "Heather Santos - default_campaign_stochastic:index_burden_Measles-PSU-Ferrari.csv.xz",
+    "Heather Santos - default_campaign_only_stochastic:index_burden_Measles-PSU-Ferrari.csv.xz",
+    "stochastic:index_burden_Measles-PSU-Ferrari.csv.xz",
+    "Heather Santos - bestcase_campaign_only_stochastic:index_burden_Measles-PSU-Ferrari.csv.xz",
+    "Heather Santos - stop_stochastic:index_burden_Measles-PSU-Ferrari.csv.xz"),
+  index = 1:8))
+
+stub <- "stochastic_burden_estimate_measles-LSHTM-Jit"
+
+hipercow::task_create_expr(resources = hres, expr =
+stoner::stone_stochastic_standardise(
+  group = "LSHTM-Jit",
+  in_path = file.path(base_in_path, "LSHTM-Jit_Measles"),
+  out_path = file.path(base_out_path, "Measles_LSHTM-Jit"),
+  scenarios = scenarios,
+  files = c(sprintf("%s-no-vaccination_Portnoy.csv.xz", stub),
+            sprintf("%s-mcv1-default_Portnoy.csv.xz", stub),
+            sprintf("%s-mcv2-default_Portnoy.csv.xz", stub),
+            sprintf("%s-mcv1-bestcase_Portnoy.csv.xz", stub),
+            sprintf("%s-mcv2-bestcase_Portnoy.csv.xz", stub),
+            sprintf("%s-campaign-default_Portnoy.csv.xz", stub),
+            sprintf("%s-campaign-only-default_Portnoy.csv.xz", stub),
+            sprintf("%s-campaign-bestcase_Portnoy.csv.xz", stub),
+            sprintf("%s-campaign-only-bestcase_Portnoy.csv.xz", stub),
+            sprintf("%s-stop_Portnoy.csv.xz", stub))))
+
+###############
+# MenA
+
+scenarios = c("mena-no-vaccination", "mena-campaign-default",
+              "mena-campaign-bestcase", "mena-routine-default",
+              "mena-routine-bestcase")
+stub <- "Andromachi Karachaliou - stochastic-burden.201910gavi-4.MenA_Cambridge-Trotter"
+
+hipercow::task_create_expr(resources = hres, expr =
+stoner::stone_stochastic_standardise(
+  group = "Cambridge-Trotter",
+  in_path = file.path(base_in_path, "Cambridge-Trotter"),
+  out_path = file.path(base_out_path, "MenA_Cambridge-Trotter"),
+  scenarios = scenarios,
+  files = c(sprintf("%s_no-vaccination_:index.csv.xz", stub),
+            sprintf("%s_campaign-default_:index.csv.xz", stub),
+            sprintf("%s_campaign-bestcase_:index.csv.xz", stub),
+            sprintf("%s_routine-default_:index.csv.xz", stub),
+            sprintf("%s_routine-bestcase_:index.csv.xz", stub)),
+  index = 1:52))
+
+stub <- "Michael Jackson - stochastic_burden_est_MenA_KPWA"
+hipercow::task_create_expr(resources = hres, expr =
+stoner::stone_stochastic_standardise(
+  group = "KPW-Jackson",
+  in_path = file.path(base_in_path, "KPW-Jackson"),
+  out_path = file.path(base_out_path, "MenA_KPW-Jackson"),
+  scenarios = scenarios,
+  files = c(sprintf("%s_both_bestcase_:index.csv.xz", stub),
+            sprintf("%s_both_default_:index.csv.xz", stub),
+            sprintf("%s_campaign_bestcase_:index.csv.xz", stub),
+            sprintf("%s_campaign_default_:index.csv.xz", stub),
+            sprintf("%s_none_default_:index.csv.xz", stub)),
+  index = 1:26))
+
+###############
+# PCV
+
+scenarios <- c("pcv-no-vaccination", "pcv-routine-default",
+               "pcv-routine-bestcase")
+
+stoner::stone_stochastic_standardise(
+  group = "JHU-Tam",
+  in_path = file.path(base_in_path, "JHU-Tam-PCV"),
+  out_path = file.path(base_out_path, "PCV_JHU-Tam"),
+  scenarios = scenarios,
+  files = c("novac:index.csv.xz", "default:index.csv.xz", "best:index.csv.xz"),
+  index = 1:14)
+
+stoner::stone_stochastic_standardise(
+  group = "LSHTM-Clark",
+  in_path = file.path(base_in_path, "LSHTM-Clark_PCV"),
+  out_path = file.path(base_out_path, "PCV_LSHTM-Clark"),
+  scenarios = scenarios,
+  files = c("VIMC_Sp_PSA_NoVax.csv.xz", "VIMC_Sp_PSA_Default.csv.xz",
+            "VIMC_Sp_PSA_Best.csv.xz"))
+
+###############
+# Rota
+
+scenarios <- c("rota-no-vaccination", "rota-routine-default",
+               "rota-routine-bestcase")
+
+stoner::stone_stochastic_standardise(
+  group = "JHU-Tam",
+  in_path = file.path(base_in_path, "JHU-Tam-Rota"),
+  out_path = file.path(base_out_path, "Rota_JHU-Tam"),
+  scenarios = scenarios,
+  files = c("novac:index.csv.xz", "default:index.csv.xz", "best:index.csv.xz"),
+  index = 1:14)
+
+stoner::stone_stochastic_standardise(
+  group = "LSHTM-Clark",
+  in_path = file.path(base_in_path, "LSHTM-Clark_Rota"),
+  out_path = file.path(base_out_path, "Rota_LSHTM-Clark"),
+  scenarios = scenarios,
+  files = c("Hira Tanvir - VIMC_Rota_PSA_NoVax.csv.xz",
+            "Hira Tanvir - VIMC_Rota_PSA_Default.csv.xz",
+            "Hira Tanvir - VIMC_Rota_PSA_Best.csv.xz"))
+
+stoner::stone_stochastic_standardise(
+  group = "Emory-Lopman",
+  in_path = file.path(base_in_path, "Emory-Lopman"),
+  out_path = file.path(base_out_path, "Rota_Emory-Lopman"),
+  scenarios = scenarios,
+  files = paste0("Molly Steele - stochastic-burden.201910gavi-4.",
+                "Rota_Emory-Lopman_:scenario.csv.xz"))
+
+###############
+# Rubella
+
+scenarios <- c(
+  "rubella-routine-no-vaccination", "rubella-campaign-default",
+  "rubella-campaign-bestcase", "rubella-rcv1-default",
+  "rubella-rcv1-bestcase", "rubella-rcv1-rcv2-default",
+  "rubella-rcv1-rcv2-bestcase", "rubella-rcv2-default",
+  "rubella-rcv2-bestcase", "rubella-stop")
+
+stub <- "Amy Winter - stochastic_burden_est-rubella"
+
+id10 <- hipercow::task_create_expr(resources = hres, expr =
+stoner::stone_stochastic_standardise(
+  group = "JHU-Lessler",
+  in_path = file.path(base_in_path, "JHU-Lessler"),
+  out_path = file.path(base_out_path, "Rubella_JHU-Lessler"),
+  scenarios = scenarios,
+  files = c(sprintf("%s-no-vaccination_:index.csv.xz", stub),
+            sprintf("%s-campaign-default_:index.csv.xz", stub),
+            sprintf("%s-campaign-bestcase_:index.csv.xz", stub),
+            sprintf("%s-rcv1-default_:index.csv.xz", stub),
+            sprintf("%s-rcv1-bestcase_:index.csv.xz", stub),
+            sprintf("%s-rcv1-rcv2-default_:index.csv.xz", stub),
+            sprintf("%s-rcv1-rcv2-bestcase_:index.csv.xz", stub),
+            sprintf("%s-rcv2-default_:index.csv.xz", stub),
+            sprintf("%s-rcv2-bestcase_:index.csv.xz", stub),
+            sprintf("%s-stop_:index.csv.xz", stub)),
+  index = 1:12))
+
+hipercow::task_create_expr(resources = hres, expr =
+stoner::stone_stochastic_standardise(
+  group = "PHE-Vynnycky",
+  in_path = file.path(base_in_path, "PHE-Vynnycky"),
+  out_path = file.path(base_out_path, "Rubella_PHE-Vynnycky"),
+  scenarios = scenarios,
+  files = "stochastic_burden_est_:scenario_country:index.csv.xz",
+  index = 1:112))
+
+
+###############
+# Typhoid
+
+scenarios = c("typhoid-no-vaccination",
+                "typhoid-campaign-default", "typhoid-routine-default")
+
+stoner::stone_stochastic_standardise(
+  group = "IVI-Kim",
+  in_path = file.path(base_in_path, "IVI-Kim-Typhoid"),
+  out_path = file.path(base_out_path, "Typhoid_IVI-Kim"),
+  scenarios = scenarios,
+  files = c("Jong-Hoon Kim - stoch_Typhoid_novacc.csv.xz",
+           "Jong-Hoon Kim - stoch_Typhoid_campaign.csv.xz",
+           "Jong-Hoon Kim - stoch_Typhoid_routine.csv.xz"))
+
+stoner::stone_stochastic_standardise(
+  group = "Yale-Pitzer",
+  in_path = file.path(base_in_path, "Yale-Pitzer"),
+  out_path = file.path(base_out_path, "Typhoid_Yale-Pitzer"),
+  scenarios = scenarios,
+  files = c("Virginia Pitzer - 2021-02-18 17.00.26 - stochastic_output_TF-Yale-Pitzer_novacc.csv.xz",
+            "Virginia Pitzer - 2021-02-18 16.58.03 - stochastic_output_TF-Yale-Pitzer_campaign.csv.xz",
+            "Virginia Pitzer - 2021-02-18 16.59.14 - stochastic_output_TF-Yale-Pitzer_camproutine.csv.xz"))
+
+###############
+# YF
+
+scenarios <- c("yf-no-vaccination", "yf-preventive-bestcase",
+               "yf-preventive-default", "yf-routine-bestcase",
+               "yf-routine-default", "yf-stop")
+
+hipercow::task_create_expr(resources = hres, expr =
+stoner::stone_stochastic_standardise(
+  group = "IC-Garske",
+  in_path = file.path(base_in_path, "IC-Garske"),
+  out_path = file.path(base_out_path, "YF_IC-Garske"),
+  scenarios = scenarios,
+  files = "stochastic-burden-estimates.201910gavi-4_YF_IC-Garske_:scenario_:index.csv.xz",
+  index = 1:200))
+
+hipercow::task_create_expr(resources = hres, expr =
+stoner::stone_stochastic_standardise(
+  group = "UND-Perkins",
+  in_path = file.path(base_in_path, "UND-Perkins"),
+  out_path = file.path(base_out_path, "YF_UND-Perkins"),
+  scenarios = scenarios,
+  files = "stochastic_burden_est_YF_UND-Perkins_:scenario_:index.csv.xz",
+  index = 1:200))
diff --git a/scripts/process_stoch_202110gavi.R b/scripts/process_stoch_202110gavi.R
new file mode 100644
index 0000000..210d016
--- /dev/null
+++ b/scripts/process_stoch_202110gavi.R
@@ -0,0 +1,386 @@
+base_in_path <- "//wpia-hn2.hpc.dide.ic.ac.uk/vimc_stochastics_dropbox/latest/202110gavi"
+base_out_path <- "//wpia-hn2.hpc.dide.ic.ac.uk/vimc_stochastics/202110gavi"
+
+vault <- vaultr::vault_client(login = "github")
+password <- vault$read("/secret/vimc/database/production/users/readonly")$password
+con <- DBI::dbConnect(RPostgres::Postgres(),
+                      dbname = "montagu",
+                      host = "montagu.vaccineimpact.org",
+                      port = 5432, password = password,
+                      user = "readonly")
+
+fetch_scenarios <- function(disease) {
+  sort(DBI::dbGetQuery(con, "
+    SELECT scenario_description FROM
+           scenario JOIN scenario_description
+           ON scenario.scenario_description = scenario_description.id
+           WHERE disease = $1 AND touchstone='202110gavi-3'", disease)$scenario_description)
+}
+
+# Let's unleash the cow
+
+setwd("Q:/testcow")
+writeLines("vimc/stoner@VIMC-9230", "pkgdepends.txt")
+hipercow::hipercow_init(driver = "dide-windows")
+hipercow::hipercow_provision()
+# Network/memory might be too much for more than a job per node.
+hres <- hipercow::hipercow_resources(cores = 1L, exclusive = TRUE)
+
+###############
+# Cholera
+
+scenarios = c("cholera-no-vaccination", "cholera-campaign-default")
+
+hipercow::task_create_expr(resources = hres, expr =
+stoner::stone_stochastic_standardise(
+  group = "IVI-Kim",
+  in_path = file.path(base_in_path, "IVI-Kim-Cholera"),
+  out_path = file.path(base_out_path, "Cholera_IVI-Kim"),
+  scenarios = scenarios,
+  files = c("Jong-Hoon Kim - stoch_Cholera_campaign_20211221T00.csv.xz",
+            "Jong-Hoon Kim - stoch_Cholera_campaign_20211222T212131.csv.xz")
+  ))
+
+stoner::stone_stochastic_standardise(
+  group = "JHU-Lee",
+  in_path = file.path(base_in_path, "JHU-Lee-Cholera"),
+  out_path = file.path(base_out_path, "Cholera_JHU-Lee"),
+  scenarios = scenarios,
+  files = c("Kaiyue Zou - no-vaccination.csv.xz",
+            "Kaiyue Zou - campaign-default.csv.xz"))
+
+###############
+# HepB
+
+scenarios <- c("hepb-no-vaccination",
+               "hepb-bd-routine-default",
+               "hepb-bd-routine-ia2030_target",
+               "hepb-bd-routine-ia2030_target-hepb-routine-ia2030_target",
+               "hepb-bd-default-hepb-routine-default",
+               "hepb-hepb-routine-default",
+               "hepb-hepb-routine-ia2030_target")
+
+hipercow::task_create_expr(resources = hres, expr =
+  stoner::stone_stochastic_standardise(
+    group = "Li",
+    in_path = file.path(base_in_path, "Li"),
+    out_path = file.path(base_out_path, "HepB_Li"),
+    scenarios = scenarios,
+    files = ":scenario:index.csv.xz",
+    index = 1:200))
+
+hipercow::task_create_expr(resources = hres, expr =
+  stoner::stone_stochastic_standardise(
+    group = "IC-Hallett",
+    in_path = file.path(base_in_path, "IC-Hallett"),
+    out_path = file.path(base_out_path, "HepB_IC-Hallett"),
+    scenarios = scenarios,
+    files = "stochastic_burden_est_HepB-IC-Hallett_:scenario_:index.csv.xz",
+    index = 1:200))
+
+###############
+# HIB
+
+scenarios <- c("hib-no-vaccination", "hib-routine-default",
+               "hib-routine-ia2030_target")
+
+stoner::stone_stochastic_standardise(
+  group = "JHU-Tam",
+  in_path = file.path(base_in_path, "JHU-Tam-Carter-Hib"),
+  out_path = file.path(base_out_path, "Hib_JHU-Tam"),
+  scenarios = scenarios,
+  files = c("hib-no-vaccination-LiST.csv.xz",
+            "hib-routine-default-LiST.csv.xz",
+            "hib-routine-ia2030_target-LiST.csv.xz"))
+
+stoner::stone_stochastic_standardise(
+  group = "LSHTM-Clark",
+  in_path = file.path(base_in_path, "LSHTM-Clark_Hib"),
+  out_path = file.path(base_out_path, "Hib_LSHTM-Clark"),
+  scenarios = scenarios,
+  files = c("Kaja Abbas - PSA_202110gavi-3_hib-no-vaccination.csv.xz",
+            "Kaja Abbas - PSA_202110gavi-3_hib-routine-default.csv.xz",
+            "Kaja Abbas - PSA_202110gavi-3_hib-routine-ia2030_target.csv.xz"))
+
+###############
+# HPV
+
+scenarios <- c("hpv-no-vaccination", "hpv-campaign-default",
+               "hpv-routine-default","hpv-campaign-ia2030_target",
+               "hpv-routine-ia2030_target")
+
+hipercow::task_create_expr(resources = hres, expr =
+stoner::stone_stochastic_standardise(
+  group = "Harvard-Sweet",
+  in_path = file.path(base_in_path, "Harvard-Sweet"),
+  out_path = file.path(base_out_path, "HPV_Harvard-Sweet"),
+  scenarios = scenarios,
+  files = c("Allison Portnoy - stochastic-burden-est.novacc_run_200.csv.xz",
+            "Allison Portnoy - stochastic-burden-est.coverage_202110gavi-3_hpv-campaign-default_run_:index.csv.xz",
+            "Allison Portnoy - stochastic-burden-est.coverage_202110gavi-3_hpv-routine-default_run_:index.csv.xz",
+            "Allison Portnoy - stochastic-burden-est.coverage_202110gavi-3_hpv-campaign-ia2030_target_run_:index.csv.xz",
+            "Allison Portnoy - stochastic-burden-est.coverage_202110gavi-3_hpv-routine-ia2030_target_run_:index.csv.xz"),
+  index = 1:200))
+
+id <- hipercow::task_create_expr(resources = hres, expr =
+stoner::stone_stochastic_standardise(
+  group = "LSHTM-Jit",
+  in_path = file.path(base_in_path, "LSHTM-Jit_HPV"),
+  out_path = file.path(base_out_path, "HPV_LSHTM-Jit"),
+  scenarios =  scenarios,
+  files =
+    c("stochastic-burden-novaccination_all_202110gavi-3_hpv-no-vaccination.csv.xz",
+      "stochastic-burden-vaccination_all_202110gavi-3_hpv-campaign-default.csv.xz",
+      "stochastic-burden-vaccination_all_202110gavi-3_hpv-routine-default.csv.xz",
+      "stochastic-burden-vaccination_all_202110gavi-3_hpv-campaign-ia2030_target.csv.xz",
+      "stochastic-burden-vaccination_all_202110gavi-3_hpv-routine-ia2030_target.csv.xz")))
+
+###############
+# JE
+
+scenarios <- c("je-routine-no-vaccination", "je-campaign-default",
+               "je-campaign-ia2030_target", "je-routine-default",
+               "je-routine-ia2030_target")
+
+hipercow::task_create_expr(resources = hres, expr =
+stoner::stone_stochastic_standardise(
+  group = "NUS-Clapham",
+  in_path = file.path(base_in_path, "NUS-Clapham-JE"),
+  out_path = file.path(base_out_path, "JE_NUS-Clapham"),
+  scenarios =  scenarios,
+  files = c(
+    "Naive_Stochastic__:index.csv.xz",
+    "Campaign_Stochastic_:index.csv.xz",
+    "Campaign_Target_Stochastic_:index.csv.xz",
+    "Routine_Stochastic_:index.csv.xz",
+    "Routine_Target_Stochastic_:index.csv.xz"),
+  index = 1:200))
+
+stub <- "Sean Moore - stochastic_burden_est_JE_UND-Moore"
+
+hipercow::task_create_expr(resources = hres, expr =
+stoner::stone_stochastic_standardise(
+  group = "UND-Moore",
+  in_path = file.path(base_in_path, "UND-Moore-JE"),
+  out_path = file.path(base_out_path, "JE_UND-Moore"),
+  scenarios =  scenarios,
+  files = sprintf("%s_:scenario.csv.xz", stub)))
+
+
+###############
+# Measles
+scenarios <- c("measles-no-vaccination",
+               "measles-campaign-default",
+               "measles-campaign-ia2030_target",
+               "measles-campaign-only-default",
+               "measles-campaign-only-ia2030_target",
+               "measles-mcv1-default",
+               "measles-mcv1-ia2030_target",
+               "measles-mcv2-default",
+               "measles-mcv2-ia2030_target")
+
+hipercow::task_create_expr(resources = hres, expr =
+stoner::stone_stochastic_standardise(
+  group = "PSU-Ferrari",
+  in_path = file.path(base_in_path, "PSU-Ferrari-Measles"),
+  out_path = file.path(base_out_path, "Measles_PSU-Ferrari"),
+  scenarios = scenarios,
+  files = "coverage_202110gavi-3_:scenario.csv.xz"))
+
+hipercow::task_create_expr(resources = hres, expr =
+stoner::stone_stochastic_standardise(
+  group = "LSHTM-Jit",
+  in_path = file.path(base_in_path, "LSHTM-Jit_Measles"),
+  out_path = file.path(base_out_path, "Measles_LSHTM-Jit"),
+  scenarios = scenarios,
+  files = c("Han Fu - stochastic_burden_estimate_measles-LSHTM-Jit-no-vaccination.csv.xz",
+            "Han Fu - stochastic_burden_estimate_measles-LSHTM-Jit-campaign-default.csv.xz",
+            "Han Fu - stochastic_burden_estimate_measles-LSHTM-Jit-campaign-ia2030_target.csv.xz",
+            "Han Fu - stochastic_burden_estimate_measles-LSHTM-Jit-campaign-only-default.csv.xz",
+            "Han Fu - stochastic_burden_estimate_measles-LSHTM-Jit-campaign-only-ia2030_target.csv.xz",
+            "Han Fu - stochastic_burden_estimate_measles-LSHTM-Jit-mcv1-default.csv.xz",
+            "Han Fu - stochastic_burden_estimate_measles-LSHTM-Jit-mcv1-ia2030_target.csv.xz",
+            "Han Fu - stochastic_burden_estimate_measles-LSHTM-Jit-mcv2-default.csv.xz",
+            "Han Fu - stochastic_burden_estimate_measles-LSHTM-Jit-mcv2-ia2030_target.csv.xz"
+            )))
+
+###############
+# MenA
+
+scenarios = c("mena-no-vaccination", "mena-campaign-default",
+              "mena-campaign-ia2030_target", "mena-routine-default",
+              "mena-routine-ia2030_target", "mena-booster-default")
+stub <- "Andromachi Karachaliou - stochastic-burden.202110gavi-2.MenA_Cambridge-Trotter"
+
+hipercow::task_create_expr(resources = hres, expr =
+stoner::stone_stochastic_standardise(
+  group = "Cambridge-Trotter",
+  in_path = file.path(base_in_path, "Cambridge-Trotter"),
+  out_path = file.path(base_out_path, "MenA_Cambridge-Trotter"),
+  scenarios = scenarios,
+  files = c(sprintf("%s_no_vaccination_:index.csv.xz", stub),
+            sprintf("%s_campaign_default_:index.csv.xz", stub),
+            sprintf("%s_campaign-ia2030_target_:index.csv.xz", stub),
+            sprintf("%s_routine_default_:index.csv.xz", stub),
+            sprintf("%s_routine-ia2030_target_:index.csv.xz", stub),
+            sprintf("%s_booster_:index.csv.xz", stub)),
+  index = 1:26))
+
+hipercow::task_create_expr(resources = hres, expr =
+stoner::stone_stochastic_standardise(
+  group = "KPW-Jackson",
+  in_path = file.path(base_in_path, "KPW-Jackson-MenA"),
+  out_path = file.path(base_out_path, "MenA_KPW-Jackson"),
+  scenarios = scenarios,
+  files = c("stochastic_burden_est_MenA_KPWA_none_default_:index.csv.xz",
+            "stochastic_burden_est_MenA_KPWA_campaign_default_:index.csv.xz",
+            "stochastic_burden_est_MenA_KPWA_campaign_ia2030_target_:index.csv.xz",
+            "stochastic_burden_est_MenA_KPWA_routine_default_:index.csv.xz",
+            "stochastic_burden_est_MenA_KPWA_routine_ia2030_target_:index.csv.xz",
+            "stochastic_burden_est_MenA_KPWA_booster_default_:index.csv.xz"),
+  index = 1:26))
+
+###############
+# PCV
+
+scenarios <- c("pcv-no-vaccination", "pcv-routine-default",
+               "pcv-routine-ia2030_target")
+
+stoner::stone_stochastic_standardise(
+  group = "JHU-Tam",
+  in_path = file.path(base_in_path, "JHU-Tam-Carter-PCV"),
+  out_path = file.path(base_out_path, "PCV_JHU-Tam"),
+  scenarios = scenarios,
+  files = ":scenario-LiST.csv.xz")
+
+stoner::stone_stochastic_standardise(
+  group = "NUS-Chen",
+  in_path = file.path(base_in_path, "LSHTM-NUS-Chen_PCV"),
+  out_path = file.path(base_out_path, "PCV_NUS-Chen"),
+  scenarios = scenarios,
+  files = "stochastic_burden_est_:scenario.csv.xz")
+
+###############
+# Rota
+
+scenarios <- c("rota-no-vaccination", "rota-routine-default",
+               "rota-routine-ia2030_target")
+
+stoner::stone_stochastic_standardise(
+  group = "JHU-Tam",
+  in_path = file.path(base_in_path, "JHU-Tam-Carter-Rota"),
+  out_path = file.path(base_out_path, "Rota_JHU-Tam"),
+  scenarios = scenarios,
+  files = ":scenario-LiST.csv.xz")
+
+stoner::stone_stochastic_standardise(
+  group = "LSHTM-Clark",
+  in_path = file.path(base_in_path, "LSHTM-Clark_Rota"),
+  out_path = file.path(base_out_path, "Rota_LSHTM-Clark"),
+  scenarios = scenarios,
+  files = "Kaja Abbas - PSA_202110gavi-3_:scenario.csv.xz")
+
+stoner::stone_stochastic_standardise(
+  group = "Emory-Lopman",
+  in_path = file.path(base_in_path, "Emory-Lopman"),
+  out_path = file.path(base_out_path, "Rota_Emory-Lopman"),
+  scenarios = scenarios,
+  files = c("Aniruddha Deshpande - stochastic_burden_est_lopman_no_vaccination_2022_01_31.csv.xz",
+            "Aniruddha Deshpande - stochastic_burden_est_lopman_routine_2022_01_31.csv.xz",
+            "Aniruddha Deshpande - stochastic_burden_est_lopman_ia2030_target_2022_01_31.csv.xz"))
+
+###############
+# Rubella
+
+scenarios <- c("rubella-routine-no-vaccination",
+  "rubella-campaign-default", "rubella-campaign-ia2030_target",
+  "rubella-rcv1-default", "rubella-rcv1-ia2030_target",
+  "rubella-rcv1-rcv2-default", "rubella-rcv1-rcv2-ia2030_target",
+  "rubella-rcv2-default", "rubella-rcv2-ia2030_target")
+
+id10 <- hipercow::task_create_expr(resources = hres, expr =
+stoner::stone_stochastic_standardise(
+  group = "JHU-Lessler",
+  in_path = file.path(base_in_path, "JHU-UGA-Winter-Rubella"),
+  out_path = file.path(base_out_path, "Rubella_JHU-Lessler"),
+  scenarios = scenarios,
+  files = "Amy Winter - stochastic_burden_est-:scenario_:index.csv.xz",
+  index = 1:11))
+
+id11<-hipercow::task_create_expr(resources = hres, expr =
+stoner::stone_stochastic_standardise(
+  group = "PHE-Vynnycky",
+  in_path = file.path(base_in_path, "PHE-Vynnycky_Rubella"),
+  out_path = file.path(base_out_path, "Rubella_PHE-Vynnycky"),
+  scenarios = scenarios,
+  files = c("VIMC_NV_RCV1RCV2Camp_country:index.csv.xz",
+            "VIMC_DF_Camp_country:index.csv.xz",
+            "VIMC_IA_Camp_country:index.csv.xz",
+            "VIMC_DF_RCV1Camp_country:index.csv.xz",
+            "VIMC_IA_RCV1Camp_country:index.csv.xz",
+            "VIMC_DF_RCV1RCV2_country:index.csv.xz",
+            "VIMC_IA_RCV1RCV2_country:index.csv.xz",
+            "VIMC_DF_RCV1RCV2Camp_country:index.csv.xz",
+            "VIMC_IA_RCV1RCV2Camp_country:index.csv.xz"),
+  index = 1:112))
+
+
+
+###############
+# Typhoid
+
+scenarios = c("typhoid-no-vaccination",
+              "typhoid-campaign-default",
+              "typhoid-campaign-ia2030_target",
+              "typhoid-routine-default",
+              "typhoid-routine-ia2030_target")
+
+hipercow::task_create_expr(resources = hres, expr =
+stoner::stone_stochastic_standardise(
+  group = "IVI-Kim",
+  in_path = file.path(base_in_path, "IVI-Kim-Typhoid"),
+  out_path = file.path(base_out_path, "Typhoid_IVI-Kim"),
+  scenarios = scenarios,
+  files = c("Jong-Hoon Kim - stoch_Typhoid_novacc_20220509T233408.csv.xz",
+            "Jong-Hoon Kim - stoch_Typhoid_campaign-default_20220510T103619.csv.xz",
+            "Jong-Hoon Kim - stoch_Typhoid_campaign-ia2030_20220510T105521.csv.xz",
+            "Jong-Hoon Kim - stoch_Typhoid_routine-default_20220510T103357.csv.xz",
+            "Jong-Hoon Kim - stoch_Typhoid_routine-ia2030_20220514.csv.xz")))
+
+hipercow::task_create_expr(resources = hres, expr =
+stoner::stone_stochastic_standardise(
+  group = "Yale-Pitzer",
+  in_path = file.path(base_in_path, "Yale-Pitzer-Typhoid"),
+  out_path = file.path(base_out_path, "Typhoid_Yale-Pitzer"),
+  scenarios = scenarios,
+  files = c("Holly Burrows - stochastic_burden_est_TF-Yale-Burrows-novacc_202110.csv.xz",
+            "Holly Burrows - stochastic_burden_est_TF-Yale-Burrows_campaign-default_202110.csv.xz",
+            "Holly Burrows - stochastic_burden_est_TF-Yale-Burrows_campaign-IA2030_202110.csv.xz",
+            "Holly Burrows - stochastic_burden_est_TF-Yale-Burrows_routine-default_202110.csv.xz",
+            "Holly Burrows - stochastic_burden_est_TF-Yale-Burrows_routine-IA2030_202110.csv.xz")))
+
+
+
+###############
+# YF
+
+scenarios <- c("yf-no-vaccination", "yf-preventive-ia2030_target",
+               "yf-preventive-default", "yf-routine-ia2030_target",
+               "yf-routine-default")
+
+stub <- "Keith Fraser - stochastic-burden-estimates.202110gavi-3_YF_IC-Garske"
+stoner::stone_stochastic_standardise(
+  group = "IC-Garske",
+  in_path = file.path(base_in_path, "IC-Garske"),
+  out_path = file.path(base_out_path, "YF_IC-Garske"),
+  scenarios = scenarios,
+  files = sprintf("%s_:scenario_:index.csv.xz", stub),
+  index = 1:200)
+
+stoner::stone_stochastic_standardise(
+  group = "UND-Perkins",
+  in_path = file.path(base_in_path, "UND-Perkins-YF"),
+  out_path = file.path(base_out_path, "YF_UND-Perkins"),
+  scenarios = scenarios,
+  files = "stochastic_burden_est_YF_UND-Perkins_:scenario_:index.csv.xz",
+  index = 1:200)
diff --git a/scripts/process_stoch_202310gavi.R b/scripts/process_stoch_202310gavi.R
index 0c4915b..b2eaaef 100644
--- a/scripts/process_stoch_202310gavi.R
+++ b/scripts/process_stoch_202310gavi.R
@@ -20,7 +20,7 @@ stoner::stone_stochastic_standardise(
 stoner::stone_stochastic_standardise(
   group = "JHU-Lee",
   in_path = file.path(base_in_path, "Cholera-JHU-Lee"),
-  out_path = file.path(base_out_path, "Cholera-JHU-Lee"),
+  out_path = file.path(base_out_path, "Cholera_JHU-Lee"),
   scenarios = scenarios,
   files = c("stochastic-burden-template.202310gavi-4.Cholera_standard_template.529.no-vaccination Christina Alam.csv.xz",
             "stochastic-burden-template.202310gavi-4.Cholera_standard_template.529.ocv1-default_one Christina Alam.csv.xz",

From 31fa3ded13547e6c52e40c7def0eec83048f535e Mon Sep 17 00:00:00 2001
From: Wes Hinsley <w.hinsley@imperial.ac.uk>
Date: Thu, 26 Mar 2026 11:24:33 +0000
Subject: [PATCH 05/29] Vignette warning

---
 vignettes/Stochastics.Rmd | 2 +-
 vignettes/Usage.Rmd       | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/vignettes/Stochastics.Rmd b/vignettes/Stochastics.Rmd
index 10bafee..9be11c9 100644
--- a/vignettes/Stochastics.Rmd
+++ b/vignettes/Stochastics.Rmd
@@ -4,7 +4,7 @@ author: "Wes Hinsley"
 date: "`r Sys.Date()`"
 output: rmarkdown::html_vignette
 vignette: >
-  %\VignetteIndexEntry{stoner}
+  %\VignetteIndexEntry{Stochastics}
   %\VignetteEngine{knitr::rmarkdown}
   %\VignetteEncoding{UTF-8}
 ---
diff --git a/vignettes/Usage.Rmd b/vignettes/Usage.Rmd
index 89d8b17..5268108 100644
--- a/vignettes/Usage.Rmd
+++ b/vignettes/Usage.Rmd
@@ -4,7 +4,7 @@ author: "Wes Hinsley"
 date: "`r Sys.Date()`"
 output: rmarkdown::html_vignette
 vignette: >
-  %\VignetteIndexEntry{stoner}
+  %\VignetteIndexEntry{Usage}
   %\VignetteEngine{knitr::rmarkdown}
   %\VignetteEncoding{UTF-8}
 ---

From 6e680a340d1b402db9ff2788d06990d5b07ae5bc Mon Sep 17 00:00:00 2001
From: Wes Hinsley <w.hinsley@imperial.ac.uk>
Date: Thu, 26 Mar 2026 11:58:44 +0000
Subject: [PATCH 06/29] Packit, and more graph control

---
 R/stochastic_graphs.R | 88 +++++++++++++++++++++++++++++++++++++++----
 1 file changed, 81 insertions(+), 7 deletions(-)

diff --git a/R/stochastic_graphs.R b/R/stochastic_graphs.R
index eee7184..1a6e2e2 100644
--- a/R/stochastic_graphs.R
+++ b/R/stochastic_graphs.R
@@ -21,23 +21,65 @@
 ##' @param by_cohort If TRUE, then age is subtracted from year to convert it to
 ##' year of birth before aggregating.
 ##' @param log If TRUE, then use a logged y-axis.
+##' @param packit_id If set, then read central burden estimates from a file
+##' within a packit on the Montagu packit server.
+##' @param packit_file Used with packit_id to specify the filename of an RDS
+##' file providing burden estimates. We expect to find scenario, year, age,
+##' country, burden_outcome and value fields in the table.
+##' @param include_quantiles Default TRUE, select whether to plot the
+##' 5% and 95% quantile lines.
+##' @param include_mean Default TRUE, select whether to plot the mean.
+##' @param include_median Default TRUE, select whether to plot the median.
+##' @param scenario2 Default NULL; if set, then the burdens from this
+##' scenario will be subtracted from those in `scenario` - ie, this plots
+##' an impact graph of applying the second scenario. For many graphs that
+##' use this, the result will be positive numbers, representing cases
+##' or deaths averted.
 
 stone_stochastic_graph <- function(base, touchstone, disease, group, country,
                                    scenario, outcome, ages = NULL,
-                                   by_cohort = FALSE, log = FALSE) {
+                                   by_cohort = FALSE, log = FALSE,
+                                   packit_id = NULL, packit_file = NULL,
+                                   include_quantiles = TRUE,
+                                   include_mean = TRUE,
+                                   include_median = TRUE,
+                                   scenario2 = NULL) {
 
   d <- prepare_graph_data(base, touchstone, disease, group, country,
                          scenario, outcome, ages, by_cohort)
+  if (is.null(ages)) {
+    age <- "all ages"
+  } else if (all(sort(ages) == min(ages):max(ages))) {
+    age <- sprintf("age %d..%d", min(ages), max(ages))
+  } else age <- "selected ages"
 
-  title <- sprintf("%s, %s, %s\n%s, %s\n", touchstone, disease, group,
+  title <- sprintf("%s, %s, %s, %s\n%s, %s\n", touchstone, disease, group, age,
                    scenario, country)
+  outcome_ylab <- outcome
+
+  if (!is.null(scenario2)) {
+    d2 <- prepare_graph_data(base, touchstone, disease, group, country,
+                            scenario2, outcome, ages, by_cohort)
+    title <- sprintf("%s, %s, %s, %s\n%s, %s\n", touchstone, disease, group, age,
+                     sprintf("Impact of %s ->\n%s", scenario, scenario2), country)
+    d <- d[order(d$year, d$run_id), ]
+    d2 <- d2[order(d2$year, d2$run_id), ]
+    d[[outcome]] <- d[[outcome]] - d2[[outcome]]
+    outcome_ylab <- paste(outcome_ylab, "averted")
+  }
 
   runs <- max(d$run_id)
   miny <- max(1, min(d[[outcome]]))
   maxy <- max(d[[outcome]])
   log <- if (log) "y" else ""
 
-  plot(ylab = outcome, xlab = if (by_cohort) "Birth Cohort" else "year",
+  if (!is.null(packit_id)) {
+    central <- prepare_central_data(packit_id, packit_file)
+
+
+  }
+  par(mar = c(5, 4, 5, 2))
+  plot(ylab = outcome_ylab, xlab = if (by_cohort) "Birth Cohort" else "year",
        x = d$year[d$run_id == 1], y = d[[outcome]][d$run_id == 1], type="l",
        col = "#b0b0b0", ylim = c(miny, maxy), main = title, log = log)
 
@@ -54,10 +96,17 @@ stone_stochastic_graph <- function(base, touchstone, disease, group, country,
       q95 = quantile(.data[[outcome]], 0.95),
       .groups = "drop"
     )
-  lines(x = avgs$year, y = avgs$mean, col = "#ff4040", lwd = 2)
-  lines(x = avgs$year, y = avgs$median, col = "#00ff00", lwd = 2)
-  lines(x = avgs$year, y = avgs$q05, col = "#202020", lwd = 2)
-  lines(x = avgs$year, y = avgs$q95, col = "#202020", lwd = 2)
+  if (include_mean) {
+    lines(x = avgs$year, y = avgs$mean, col = "#ff4040", lwd = 2)
+  }
+  if (include_median) {
+    lines(x = avgs$year, y = avgs$median, col = "#00ff00", lwd = 2)
+  }
+  if (include_quantiles) {
+    lines(x = avgs$year, y = avgs$q05, col = "#202020", lwd = 2)
+    lines(x = avgs$year, y = avgs$q95, col = "#202020", lwd = 2)
+  }
+  recordPlot()
 }
 
 
@@ -82,3 +131,28 @@ prepare_graph_data <- function(base, touchstone, disease, group, country,
   d
 }
 
+prepare_central_data <- function(packit_id, packit_file,
+  country, scenario, outcome, ages, by_cohort) {
+
+  central <- readRDS(fetch_packit(packit_id, packit_file))
+  central <- central[central$country == country, ]
+  central <- central[central$scenario == scenario, ]
+  central <- central[central$burden_outcome == outcome, ]
+  if (!is.null(ages)) {
+    central <- central[d$age %in% ages, ]
+  }
+  if (by_cohort) {
+    central$year <- central$year - d$age
+  }
+  names(central)[names(central) == "value"] <- outcome
+  central <- central %>% group_by(.data$year) %>%
+    summarise(
+      !!outcome := sum(.data[[outcome]], na.rm = TRUE),
+      .groups = "drop")
+  central$run_id <- 0
+  central
+}
+
+run_app <- function() {
+  shiny::runApp(system.file("app", package = "stoner"))
+}

From 65828b6e6a959d6161e2558ed9b5636755ca63ba Mon Sep 17 00:00:00 2001
From: Wes Hinsley <w.hinsley@imperial.ac.uk>
Date: Thu, 26 Mar 2026 12:00:54 +0000
Subject: [PATCH 07/29] Redocument

---
 NAMESPACE                           |  3 ++
 R/packit.R                          | 78 +++++++++++++++++++++++++++++
 man/fetch_packit.Rd                 | 29 +++++++++++
 man/stone_stochastic_graph.Rd       | 28 ++++++++++-
 man/stone_stochastic_standardise.Rd | 24 ++++++++-
 5 files changed, 160 insertions(+), 2 deletions(-)
 create mode 100644 R/packit.R
 create mode 100644 man/fetch_packit.Rd

diff --git a/NAMESPACE b/NAMESPACE
index d71f8c5..e07131a 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -1,5 +1,6 @@
 # Generated by roxygen2: do not edit by hand
 
+export(fetch_packit)
 export(stone_dump)
 export(stone_extract)
 export(stone_load)
@@ -16,6 +17,7 @@ export(stoner_calculate_dalys)
 export(stoner_dalys_for_db)
 import(arrow)
 import(dplyr)
+import(httr2)
 import(readr)
 importFrom(arrow,write_parquet)
 importFrom(data.table,as.data.table)
@@ -27,6 +29,7 @@ importFrom(stats,median)
 importFrom(stats,quantile)
 importFrom(testthat,expect_equal)
 importFrom(testthat,expect_true)
+importFrom(tools,file_ext)
 importFrom(utils,capture.output)
 importFrom(utils,read.csv)
 importFrom(utils,write.csv)
diff --git a/R/packit.R b/R/packit.R
new file mode 100644
index 0000000..03177cd
--- /dev/null
+++ b/R/packit.R
@@ -0,0 +1,78 @@
+##' Download a file from a packit, for example an artefact containing
+##' central burden estimates, which we could then plot on top of
+##' stochastics. But more generally, we can use this for fetching any
+##' file from packet (ie, from the Montagu Reporting Portal).
+##'
+##' @export
+##' @title Fetch packit
+##' @import httr2
+##' @importFrom tools file_ext
+##' @param packit_id The id of the packit containing the artefact.
+##' @param filename The filename of the file within the packit.
+##' @param server By default, the URL to the packit API on Montagu,
+##' but this can be set to other packit API's if we want.
+##' @returns The filename of the temporary file which has been downloaded.
+
+fetch_packit <- function(packet_id, filename,
+  server = "https://montagu.vaccineimpact.org/packit/api/") {
+
+  # First we have to create a client, and a flow...
+
+  client <- httr2::oauth_client(
+    id = "orderly",
+    token_url = sprintf("%s/deviceAuth/token", server),
+    name = "orderly"
+  )
+
+  # You will be asked to type a code at this point...
+
+  flow <- httr2::oauth_flow_device(
+    client = client,
+    auth_url = sprintf("%s/deviceAuth", server),
+    pkce = FALSE,
+    scope = NULL,
+    open_browser = FALSE,
+    auth_params = list(),
+    token_params = list()
+  )
+
+  # Now we have a flow$access_token we can use, to get
+  # a one-time token to download the file.
+
+  req <- httr2::request(sprintf("%s/packets/%s/files/token", server, packet_id))
+  req <- req |>
+    httr2::req_headers("Accept" = "application/json",
+                       "Content-Type" = "application/json",
+                       "Authorization" = paste("Bearer", flow$access_token)) |>
+         httr2::req_method("POST") |>
+         httr2::req_body_json(list(paths = list(filename)))
+
+  x <- httr2::req_perform(req)
+  ott <- httr2::resp_body_json(x)$id
+
+  # This is the URL to the file we want, including the token.
+
+  url <- sprintf("%s/packets/%s/file?path=%s&token=%s&filename=%s&inline=false",
+                 server, packet_id, filename, ott, filename)
+
+  req <- httr2::request(url)
+  req <- req |>
+    httr2::req_headers("Authorization" = paste("Bearer", flow$access_token))
+
+
+  # And finally, we download the file in chunks.
+
+  out <- tempfile(fileext = tools::file_ext(filename))
+  outfile <- file(out, open = "wb")
+  con <- httr2::req_perform_connection(req, blocking = TRUE)
+  while (!httr2::resp_stream_is_complete(con)) {
+    chunk <- httr2::resp_stream_raw(con, kb = 32)
+    if (length(chunk) == 0) break
+    writeBin(chunk, outfile)
+  }
+  close(outfile)
+  out
+}
+
+
+
diff --git a/man/fetch_packit.Rd b/man/fetch_packit.Rd
new file mode 100644
index 0000000..3326eca
--- /dev/null
+++ b/man/fetch_packit.Rd
@@ -0,0 +1,29 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/packit.R
+\name{fetch_packit}
+\alias{fetch_packit}
+\title{Fetch packit}
+\usage{
+fetch_packit(
+  packet_id,
+  filename,
+  server = "https://montagu.vaccineimpact.org/packit/api/"
+)
+}
+\arguments{
+\item{filename}{The filename of the file within the packit.}
+
+\item{server}{By default, the URL to the packit API on Montagu,
+but this can be set to other packit API's if we want.}
+
+\item{packit_id}{The id of the packit containing the artefact.}
+}
+\value{
+The filename of the temporary file which has been downloaded.
+}
+\description{
+Download a file from a packit, for example an artefact containing
+central burden estimates, which we could then plot on top of
+stochastics. But more generally, we can use this for fetching any
+file from packet (ie, from the Montagu Reporting Portal).
+}
diff --git a/man/stone_stochastic_graph.Rd b/man/stone_stochastic_graph.Rd
index 9fa15c9..5e40693 100644
--- a/man/stone_stochastic_graph.Rd
+++ b/man/stone_stochastic_graph.Rd
@@ -14,7 +14,13 @@ stone_stochastic_graph(
   outcome,
   ages = NULL,
   by_cohort = FALSE,
-  log = FALSE
+  log = FALSE,
+  packit_id = NULL,
+  packit_file = NULL,
+  include_quantiles = TRUE,
+  include_mean = TRUE,
+  include_median = TRUE,
+  scenario2 = NULL
 )
 }
 \arguments{
@@ -40,6 +46,26 @@ if left as NULL, then all ages are used and aggregated.}
 year of birth before aggregating.}
 
 \item{log}{If TRUE, then use a logged y-axis.}
+
+\item{packit_id}{If set, then read central burden estimates from a file
+within a packit on the Montagu packit server.}
+
+\item{packit_file}{Used with packit_id to specify the filename of an RDS
+file providing burden estimates. We expect to find scenario, year, age,
+country, burden_outcome and value fields in the table.}
+
+\item{include_quantiles}{Default TRUE, select whether to plot the
+5\% and 95\% quantile lines.}
+
+\item{include_mean}{Default TRUE, select whether to plot the mean.}
+
+\item{include_median}{Default TRUE, select whether to plot the median.}
+
+\item{scenario2}{Default NULL; if set, then the burdens from this
+scenario will be subtracted from those in \code{scenario} - ie, this plots
+an impact graph of applying the second scenario. For many graphs that
+use this, the result will be positive numbers, representing cases
+or deaths averted.}
 }
 \description{
 Draw a stochastic plot showing all the different runs, with the mean,
diff --git a/man/stone_stochastic_standardise.Rd b/man/stone_stochastic_standardise.Rd
index 343275f..dcf42ae 100644
--- a/man/stone_stochastic_standardise.Rd
+++ b/man/stone_stochastic_standardise.Rd
@@ -12,7 +12,12 @@ stone_stochastic_standardise(
   files,
   index = 1,
   rubella_fix = TRUE,
-  missing_run_id_fix = TRUE
+  hepb2019_fix = TRUE,
+  hib2019_fix = TRUE,
+  missing_run_id_fix = TRUE,
+  allow_missing_yll = TRUE,
+  allow_missing_dalys = TRUE,
+  allow_missing_indexes = FALSE
 )
 }
 \arguments{
@@ -47,10 +52,27 @@ these to the simpler names. Processing Rubella stochastic files without
 this set to TRUE will fail - so while we should always do this, keeping
 the parameter makes it more clear in the code what we're doing and why.}
 
+\item{hepb2019_fix}{In 2019, HepB deaths and cases were subdivided into a
+number of different causes. This flag combines those into the single
+appropriate burden outcome.}
+
+\item{hib2019_fix}{In 2019, HepB deaths and cases were subdivided into a
+number of different causes. This flag combines those into the single
+appropriate burden outcome.}
+
 \item{missing_run_id_fix}{Some groups in the past have omitted run_id
 from the files, but included them in the filenames. This fix inserts
 that into the files if the index parameter indicates we have 200 runs to
 process.}
+
+\item{allow_missing_yll}{yll was introduced in 2023? This flag allows
+it to be missing for processing older stochastics.}
+
+\item{allow_missing_indexes}{In some early runs, different groups
+provided different numbers of files for each scenario, because some
+countries did not implement particular coverage campaigns. This
+flag needs to be TRUE for those groups, but the default is FALSE,
+since it's rare, and we generally want errors for missing files.}
 }
 \description{
 Convert a modelling group's stochastic files into an intermediate

From fbb32ce6ed44920fffea46731756cd6b476080fe Mon Sep 17 00:00:00 2001
From: Wes Hinsley <w.hinsley@imperial.ac.uk>
Date: Thu, 26 Mar 2026 12:01:19 +0000
Subject: [PATCH 08/29] Stochastic explorer!

---
 inst/app/app.R   | 427 +++++++++++++++++++++++++++++++++++++++++++++++
 inst/app/utils.R | 270 ++++++++++++++++++++++++++++++
 2 files changed, 697 insertions(+)
 create mode 100644 inst/app/app.R
 create mode 100644 inst/app/utils.R

diff --git a/inst/app/app.R b/inst/app/app.R
new file mode 100644
index 0000000..1bc368b
--- /dev/null
+++ b/inst/app/app.R
@@ -0,0 +1,427 @@
+library(shiny)
+library(shinyjs)
+library(arrow)
+
+data_dir <- "//wpia-hn2.hpc.dide.ic.ac.uk/vimc_stochastics"
+
+source("utils.R")
+
+app_ui <- function() {
+  useShinyjs()
+  fluidPage(
+    tags$style(HTML("
+      .selectize-input {
+        white-space: nowrap;
+        overflow: hidden;
+        text-overflow: ellipsis;
+      }
+
+      .shiny-options-group {
+        display: grid;
+        grid-template-columns: repeat(2, 1fr);  /* 2 per row */
+        gap: 4px 12px;  /* row gap, column gap */
+      }
+
+      .shiny-options-group .checkbox label {
+        display: flex;
+        align-items: centre;
+        margin: 0;
+        padding-left: 0;
+        gap: 6px;
+      }
+
+      .control-label {
+        margin-bottom: 2px;
+        font-size: 12px;
+      }
+
+      .selectize-input {
+        min-height: 30px;
+        padding: 4px 8px;
+        font-size: 13px;
+      }
+
+      .form-group {
+        margin-bottom: 8px;
+      }
+
+      .shiny-plot-output {
+        padding-top: 5px;
+      }
+    ")),
+    titlePanel("Stochastic Explorer"),
+    tabsetPanel(
+      tabPanel("Burden",
+        sidebarLayout(
+          sidebarPanel(
+            selectInput("b_touchstone", "Touchstone:",
+                        choices = get_touchstones(), selected = NULL),
+            selectInput("b_disease", "Disease:",
+                        choices = character(0)),
+            selectInput("b_group", "Group:",
+                        choices = character(0)),
+            selectInput("b_scenario", "Scenario:",
+                        choices = character(0)),
+            selectInput("b_country", "Country:",
+                        choices = character(0)),
+            selectInput("b_outcome", "Y-Axis:",
+                        choices = character(0)),
+            radioButtons("b_year", "X-Axis:", inline = TRUE,
+                         choices = c("Calendar", "Cohort")),
+            radioButtons("b_ages", "Age:", inline = TRUE,
+                         choices = c("All", "Under 5")),
+            checkboxGroupInput(
+              "b_options", "Plot Options:",
+              choices = c("Quantiles", "Median", "Mean", "Log-Y"),
+              selected = c("Quantiles", "Median", "Mean", "Log-Y"),
+              inline = TRUE),
+            actionButton("b_plot_btn", "Plot")
+          ), mainPanel(
+            plotOutput("b_main_plot")
+          )
+        )
+      ),
+      tabPanel("Burden/TS",
+        sidebarLayout(
+          sidebarPanel(
+            selectInput("bts_touchstone1", "Touchstone 1:",
+              choices = get_touchstones(), selected = NULL),
+            selectInput("bts_touchstone2", "Touchstone 2:",
+              choices = get_touchstones(), selected = NULL),
+            selectInput("bts_disease", "Disease:",
+              choices = character(0)),
+            selectInput("bts_group", "Group:",
+              choices = character(0)),
+            selectInput("bts_scenario", "Scenario:",
+              choices = character(0)),
+            selectInput("bts_country", "Country:",
+              choices = character(0)),
+            selectInput("bts_outcome", "Y-Axis:",
+              choices = character(0)),
+            radioButtons("bts_year", "X-Axis:", inline = TRUE,
+              choices = c("Calendar", "Cohort")),
+            radioButtons("bts_ages", "Age:", inline = TRUE,
+              choices = c("All", "Under 5")),
+            checkboxGroupInput("bts_options", "Plot Options:",
+              choices = c("Quantiles", "Median", "Mean", "Log-Y"),
+                       selected = c("Quantiles", "Median", "Mean", "Log-Y"),
+                       inline = TRUE),
+            actionButton("bts_plot_btn", "Plot")
+          ), mainPanel(
+            plotOutput("bts_main_plot1"),
+            plotOutput("bts_main_plot2")
+          )
+        )
+      ),
+      tabPanel("Burden/MG",
+        sidebarLayout(
+          sidebarPanel(
+            selectInput("bmg_touchstone", "Touchstone:",
+                        choices = get_touchstones(), selected = NULL),
+            selectInput("bmg_disease", "Disease:",
+                        choices = character(0)),
+            selectInput("bmg_group1", "Group 1:",
+                        choices = character(0)),
+            selectInput("bmg_group2", "Group 2:",
+                        choices = character(0)),
+            selectInput("bmg_scenario", "Scenario:",
+                        choices = character(0)),
+            selectInput("bmg_country", "Country:",
+                        choices = character(0)),
+            selectInput("bmg_outcome", "Y-Axis:",
+                        choices = character(0)),
+            radioButtons("bmg_year", "X-Axis:", inline = TRUE,
+                         choices = c("Calendar", "Cohort")),
+            radioButtons("bmg_ages", "Age:", inline = TRUE,
+                         choices = c("All", "Under 5")),
+            checkboxGroupInput("bmg_options", "Plot Options:",
+                         choices = c("Quantiles", "Median", "Mean", "Log-Y"),
+                         selected = c("Quantiles", "Median", "Mean", "Log-Y"),
+                         inline = TRUE),
+            actionButton("bmg_plot_btn", "Plot")
+          ), mainPanel(
+            plotOutput("bmg_main_plot1"),
+            plotOutput("bmg_main_plot2")
+          )
+        )
+      ),
+      tabPanel("Impact",
+        sidebarLayout(
+          sidebarPanel(
+            selectInput("i_touchstone", "Touchstone:",
+                       choices = get_touchstones(), selected = NULL),
+            selectInput("i_disease", "Disease:",
+                       choices = character(0)),
+            selectInput("i_group", "Group:",
+                       choices = character(0)),
+            selectInput("i_scenario1", "Base Scenario:",
+                       choices = character(0)),
+            selectInput("i_scenario2", "Compare Scenario:",
+                       choices = character(0)),
+            selectInput("i_country", "Country:",
+                       choices = character(0)),
+            selectInput("i_outcome", "Y-Axis:",
+                       choices = character(0)),
+            radioButtons("i_year", "X-Axis:", inline = TRUE,
+                        choices = c("Calendar", "Cohort")),
+            radioButtons("i_ages", "Age:", inline = TRUE,
+                        choices = c("All", "Under 5")),
+            checkboxGroupInput(
+                 "i_options", "Plot Options:",
+                 choices = c("Quantiles", "Median", "Mean", "Log-Y"),
+                 selected = c("Quantiles", "Median", "Mean", "Log-Y"),
+                 inline = TRUE),
+            actionButton("i_plot_btn", "Plot")
+          ), mainPanel(
+           plotOutput("i_main_plot")
+          )
+        )
+      )
+    )
+  )
+}
+
+app_server <- function(input, output, session) {
+
+  # Events for the basic burden tab
+
+  observeEvent(input$b_touchstone, {
+    update_touchstone(session, input$b_touchstone, input) })
+  observeEvent(input$b_disease, {
+    update_disease(session, input$b_touchstone, input$b_disease, input) })
+  observeEvent(input$b_group, {
+    update_group(session, input$b_touchstone, input$b_disease, input$b_group, input) })
+  observeEvent(input$b_scenario, {
+    update_scenario(session, input$b_touchstone, input$b_disease, input$b_group,
+                    input$b_scenario, input) })
+  observeEvent(input$b_country, {
+    update_country(session, input$b_touchstone, input$b_disease, input$b_group,
+                   input$b_scenario, input$b_country, input) })
+
+
+  plot_obj_b <- eventReactive(input$b_plot_btn, {
+    ages <- NULL
+    if (input$b_ages == "Under 5") {
+      ages <- 0:4
+    }
+    stone_stochastic_graph(
+      base = data_dir,
+      touchstone = input$b_touchstone,
+      disease = input$b_disease,
+      group = input$b_group,
+      country = input$b_country,
+      scenario = input$b_scenario,
+      outcome = input$b_outcome,
+      by_cohort = input$b_year == "Cohort",
+      ages = ages,
+      log = "Log-Y" %in% input$b_options,
+      include_median = "Median" %in% input$b_options,
+      include_quantiles = "Quantiles" %in% input$b_options,
+      include_mean = "Mean" %in% input$b_options,
+    )
+  })
+
+  output$b_main_plot <- renderPlot({
+    plot_obj_b()
+  })
+
+  # Events for the impact tab
+
+  observeEvent(input$i_touchstone, {
+    update_touchstone(session, input$i_touchstone, input, "i") })
+  observeEvent(input$i_disease, {
+    update_disease(session, input$i_touchstone, input$i_disease, input, "i") })
+  observeEvent(input$i_group, {
+    update_group(session, input$i_touchstone, input$i_disease,
+                 input$i_group, input, "i", "1")
+    update_group(session, input$i_touchstone, input$i_disease,
+                 input$i_group, input, "i", "2")
+  })
+  observeEvent(input$i_scenario1, {
+    update_scenario(session, input$i_touchstone, input$i_disease, input$i_group,
+                    input$i_scenario1, input, "i",
+                    extra_scenario = input$i_scenario2) })
+  observeEvent(input$i_scenario2, {
+    update_scenario(session, input$i_touchstone, input$i_disease, input$i_group,
+                    input$i_scenario2, input, "i",
+                    extra_scenario = input$i_scenario2) })
+  observeEvent(input$i_country, {
+    update_country(session, input$i_touchstone, input$i_disease, input$i_group,
+                   input$i_scenario1, input$i_country, input, "i") })
+
+  plot_obj_i <- eventReactive(input$i_plot_btn, {
+
+    if (input$i_scenario1 == input$i_scenario2) {
+      showModal(modalDialog(
+        title = "Warning",
+        "Same scenario selected - impact will be zero.",
+        easyClose = TRUE,
+        footer = modalButton("OK")
+      ))
+      return(NULL)
+    }
+
+    ages <- NULL
+    if (input$i_ages == "Under 5") {
+      ages <- 0:4
+    }
+    stone_stochastic_graph(
+      base = data_dir,
+      touchstone = input$i_touchstone,
+      disease = input$i_disease,
+      group = input$i_group,
+      country = input$i_country,
+      scenario = input$i_scenario1,
+      scenario2 = input$i_scenario2,
+      outcome = input$i_outcome,
+      by_cohort = input$i_year == "Cohort",
+      ages = ages,
+      log = "Log-Y" %in% input$i_options,
+      include_median = "Median" %in% input$i_options,
+      include_quantiles = "Quantiles" %in% input$i_options,
+      include_mean = "Mean" %in% input$i_options,
+    )
+  })
+
+  output$i_main_plot <- renderPlot({
+    plot_obj_i()
+  })
+
+  # Events for comparing burdens on different touchstones
+
+  observeEvent(input$bts_touchstone1, {
+    update_touchstone_ts(session, input$bts_touchstone1, input$bts_touchstone2,
+                         input, "bts") })
+  observeEvent(input$bts_touchstone2, {
+    update_touchstone_ts(session, input$bts_touchstone1, input$bts_touchstone2,
+                         input, "bts") })
+  observeEvent(input$bts_disease, {
+    update_disease_ts(session, input$bts_touchstone1, input$bts_touchstone2,
+                      input$bts_disease, input, "bts") })
+  observeEvent(input$bts_group, {
+    update_group_ts(session, input$bts_touchstone1, input$bts_touchstone2,
+                    input$bts_disease, input$bts_group, input, "bts")
+  })
+  observeEvent(input$bts_scenario, {
+    update_scenario_ts(session, input$bts_touchstone1, input$bts_touchstone2,
+                       input$bts_disease, input$bts_group, input$bts_scenario,
+                       input, "bts") })
+  observeEvent(input$bts_country, {
+    update_country_ts(session, input$bts_touchstone1, input$bts_touchstone2,
+                      input$bts_disease, input$bts_group, input$bts_scenario,
+                      input$bts_country, input, "bts") })
+
+  plot_obj_bts <- eventReactive(input$bts_plot_btn, {
+    if (input$bts_touchstone1 == input$bts_touchstone2) {
+      showModal(modalDialog(
+        title = "Warning",
+        "Same touchstones selected.",
+        easyClose = TRUE,
+        footer = modalButton("OK")
+      ))
+      return(NULL)
+    }
+    ages <- NULL
+    if (input$i_ages == "Under 5") {
+      ages <- 0:4
+    }
+    plots <- list()
+    for (ts in 1:2) {
+      plots[[ts]] <-
+        stone_stochastic_graph(
+          base = data_dir,
+          touchstone = if (ts == 1) input$bts_touchstone1 else input$bts_touchstone2,
+          disease = input$bts_disease,
+          group = input$bts_group,
+          country = input$bts_country,
+          scenario = input$bts_scenario,
+          outcome = input$bts_outcome,
+          by_cohort = input$bts_year == "Cohort",
+          ages = ages,
+          log = "Log-Y" %in% input$bts_options,
+          include_median = "Median" %in% input$bts_options,
+          include_quantiles = "Quantiles" %in% input$bts_options,
+          include_mean = "Mean" %in% input$bts_options,
+        )
+    }
+    list(plot1 = plots[[1]], plot2 = plots[[2]])
+  })
+
+  output$bts_main_plot1 <- renderPlot({
+    replayPlot(plot_obj_bts()$plot1)
+  })
+
+  output$bts_main_plot2 <- renderPlot({
+    replayPlot(plot_obj_bts()$plot2)
+  })
+
+  # Compare burden estimates for different groups
+
+  observeEvent(input$bmg_touchstone, {
+    update_touchstone(session, input$bmg_touchstone, input, "bmg") })
+  observeEvent(input$bmg_disease, {
+    update_disease_mg(session, input$bmg_touchstone, input$bmg_disease,
+                   input, "bmg") })
+  observeEvent(input$bmg_group1, {
+    update_group_mg(session, input$bmg_touchstone, input$bmg_disease,
+                    input$bmg_group1, input$bmg_group2, input, "bmg")
+  })
+  observeEvent(input$bmg_group2, {
+    update_group_mg(session, input$bmg_touchstone, input$bmg_disease,
+                    input$bmg_group1, input$bmg_group2, input, "bmg")
+  })
+  observeEvent(input$bmg_scenario, {
+    update_scenario_mg(session, input$bmg_touchstone, input$bmg_disease,
+                       input$bmg_group1, input$bmg_group2, input$bmg_scenario,
+                       input, "bmg") })
+  observeEvent(input$bmg_country, {
+    update_country_mg(session, input$bmg_touchstone, input$bmg_disease,
+                      input$bmg_group1, input$bmg_group2, input$bmg_scenario,
+                      input$bmg_country, input, "bmg") })
+
+  plot_obj_bmg <- eventReactive(input$bmg_plot_btn, {
+    if (input$bmg_group1 == input$bmg_group2) {
+      showModal(modalDialog(
+        title = "Warning",
+        "Same groups selected.",
+        easyClose = TRUE,
+        footer = modalButton("OK")
+      ))
+      return(NULL)
+    }
+    ages <- NULL
+    if (input$i_ages == "Under 5") {
+      ages <- 0:4
+    }
+    plots <- list()
+    for (ts in 1:2) {
+      plots[[ts]] <-
+        stone_stochastic_graph(
+          base = data_dir,
+          touchstone = input$bmg_touchstone,
+          disease = input$bmg_disease,
+          group = if (ts == 1) input$bmg_group1 else input$bmg_group2,
+          country = input$bmg_country,
+          scenario = input$bmg_scenario,
+          outcome = input$bmg_outcome,
+          by_cohort = input$bmg_year == "Cohort",
+          ages = ages,
+          log = "Log-Y" %in% input$bmg_options,
+          include_median = "Median" %in% input$bmg_options,
+          include_quantiles = "Quantiles" %in% input$bmg_options,
+          include_mean = "Mean" %in% input$bmg_options,
+        )
+    }
+    list(plot1 = plots[[1]], plot2 = plots[[2]])
+  })
+
+  output$bmg_main_plot1 <- renderPlot({
+    replayPlot(plot_obj_bmg()$plot1)
+  })
+
+  output$bmg_main_plot2 <- renderPlot({
+    replayPlot(plot_obj_bmg()$plot2)
+  })
+}
+
+shinyApp(app_ui, app_server)
diff --git a/inst/app/utils.R b/inst/app/utils.R
new file mode 100644
index 0000000..f68393f
--- /dev/null
+++ b/inst/app/utils.R
@@ -0,0 +1,270 @@
+# Helper functions to dig information from the
+# shared file system.
+
+# Functions ending in _ts will find things common
+# between two touchstones, for comparisons.
+
+# Functions ending in _mg will find things common
+# between to modelling groups, for comparisons.
+
+data_dir <- "//wpia-hn2.hpc.dide.ic.ac.uk/vimc_stochastics"
+
+
+get_touchstones <- function() {
+  basename(list.dirs(data_dir, recursive = FALSE))
+}
+
+get_diseases <- function(touchstone) {
+  path <- file.path(data_dir, touchstone)
+  dirs <- basename(list.dirs(path, recursive = FALSE))
+  unique(unlist(lapply(strsplit(dirs, "_"), `[[`, 1)))
+}
+
+get_diseases_ts <- function(touchstone1, touchstone2) {
+  res <- get_diseases(touchstone1)
+  res[res %in% get_diseases(touchstone2)]
+}
+
+get_groups <- function(touchstone, disease) {
+  path <- file.path(data_dir, touchstone)
+  dirs <- basename(list.dirs(path, recursive = FALSE))
+  dirs <- dirs[substr(dirs, 1, nchar(disease) + 1) == paste0(disease, "_")]
+  unique(substring(dirs, nchar(disease) + 2))
+}
+
+get_groups_ts <- function(touchstone1, touchstone2, disease) {
+  res <- get_groups(touchstone1, disease)
+  res[res %in% get_groups(touchstone2, disease)]
+}
+
+get_scenarios <- function(touchstone, disease, group) {
+  path <- file.path(data_dir, touchstone, paste(disease, group, sep = "_"))
+  files <- basename(list.files(path, recursive = FALSE))
+  files <- gsub(".pq", "", files)
+  files <- substring(files, nchar(group) + 2)
+  ends <- unlist(lapply(gregexpr("_", files), `[[`, 1)) - 1
+  unique(substring(files, 1, ends))
+}
+
+get_scenarios_ts <- function(touchstone1, touchstone2, disease, group) {
+  res <- get_scenarios(touchstone1, disease, group)
+  res[res %in% get_scenarios(touchstone2, disease, group)]
+}
+
+get_scenarios_mg <- function(touchstone, disease, group1, group2) {
+  res <- get_scenarios(touchstone, disease, group1)
+  res[res %in% get_scenarios(touchstone, disease, group2)]
+}
+
+get_countries <- function(touchstone, disease, group, scenario, extra_scenario = NULL) {
+  scenarios <- c(scenario, extra_scenario)
+  final <- NULL
+  for (scenario in scenarios) {
+    path <- file.path(data_dir, touchstone, paste(disease, group, sep = "_"))
+    files <- basename(list.files(path, recursive = FALSE))
+    files <- gsub(".pq", "", files)
+    files <- substring(files, nchar(group) + 2)
+    files <- files[substring(files, 1, nchar(scenario) + 1) == paste0(scenario, "_")]
+    ends <- unlist(lapply(gregexpr("_", files), `[[`, 1)) + 1
+    if (is.null(final)) final <- unique(substring(files, ends))
+    else final <- final[final %in% unique(substring(files, ends))]
+  }
+  final
+}
+
+get_countries_ts <- function(touchstone1, touchstone2, disease, group, scenario) {
+  res <- get_countries(touchstone1, disease, group, scenario)
+  res[res %in% get_countries(touchstone2, disease, group, scenario)]
+}
+
+get_countries_mg <- function(touchstone, disease, group1, group2, scenario) {
+  res <- get_countries(touchstone, disease, group1, scenario)
+  res[res %in% get_countries(touchstone, disease, group2, scenario)]
+}
+
+get_outcomes <- function(touchstone, disease, group, scenario, country) {
+  path <- file.path(data_dir, touchstone, paste(disease, group, sep = "_"))
+  thefile <- sprintf("%s_%s_%s.pq", group, scenario, country)
+  tbl <- arrow::read_parquet(file.path(path, thefile), col_select = NULL)
+  cols <- names(tbl)
+  sort(cols[!cols %in% c("disease", "run_id", "year", "age",
+                    "country", "cohort_size")])
+}
+
+get_outcomes_ts <- function(touchstone1, touchstone2, disease, group,
+                            scenario, country) {
+  res <- get_outcomes(touchstone1, disease, group, scenario, country)
+  res[res %in% get_outcomes(touchstone2, disease, group, scenario, country)]
+}
+
+get_outcomes_mg <- function(touchstone, disease, group1, group2, scenario,
+                            country) {
+  res <- get_outcomes(touchstone, disease, group1, scenario, country)
+  res[res %in% get_outcomes(touchstone, disease, group2, scenario, country)]
+}
+
+
+# GUI helpers.
+
+# Update a dropdown. Keep the previously selected value if possible,
+# otherwise select the first one.
+
+update_dropdown_keep <- function(session, id, choices, selected) {
+  if ((is.null(selected)) ||
+      (!selected %in% choices)) selected <- choices[[1]]
+  updateSelectInput(session, id, choices = choices, selected = selected)
+  selected
+}
+
+# The next events recurse through the dropdowns updating available
+# options. In a messy attempt to avoid duplication, the prefix is
+# "b" for the simple burden panel, "i" for impact,
+# "bts" for multi-touchstone burden, "its" for multi-touchstone impact
+# "bmg" for multi-group burden, and "img" for multi-group impact.
+
+update_touchstone <- function(session, touchstone, input, prefix = "b") {
+  d <- paste(prefix, "disease", sep = "_")
+  diseases <- get_diseases(touchstone)
+  disease <- update_dropdown_keep(session, d, diseases, input[[d]])
+  update_disease(session, touchstone, disease, input, prefix)
+}
+
+update_touchstone_ts <- function(session, touchstone1, touchstone2,
+                                 input, prefix = "bts") {
+  d <- paste(prefix, "disease", sep = "_")
+  diseases <- get_diseases_ts(touchstone1, touchstone2)
+  disease <- update_dropdown_keep(session, d, diseases, input[[d]])
+  update_disease_ts(session, touchstone1, touchstone2, disease, input, prefix)
+}
+
+update_disease <- function(session, touchstone, disease, input, prefix = "b") {
+  if (is.null(disease) || (disease == "")) return()
+  g <- paste(prefix, "group", sep = "_")
+  groups <- get_groups(touchstone, disease)
+  group <- update_dropdown_keep(session, g, groups, input[[g]])
+  update_group(session, touchstone, disease, group, input, prefix)
+}
+
+update_disease_ts <- function(session, touchstone1, touchstone2, disease,
+                              input, prefix = "bts") {
+  if (is.null(disease) || (disease == "")) return()
+  g <- paste(prefix, "group", sep = "_")
+  groups <- get_groups_ts(touchstone1, touchstone2, disease)
+  group <- update_dropdown_keep(session, g, groups, input[[g]])
+  update_group_ts(session, touchstone1, touchstone2, disease, group,
+                  input, prefix)
+}
+
+update_disease_mg <- function(session, touchstone, disease, input,
+                              prefix = "bmg") {
+  if (is.null(disease) || (disease == "")) return()
+  g1 <- paste(prefix, "group1", sep = "_")
+  g2 <- paste(prefix, "group2", sep = "_")
+  groups <- get_groups(touchstone, disease)
+  group1 <- update_dropdown_keep(session, g1, groups, input[[g1]])
+  group2 <- update_dropdown_keep(session, g2, groups, input[[g2]])
+  update_group_mg(session, touchstone, disease, group1, group2,
+                  input, prefix)
+}
+
+update_group <- function(session, touchstone, disease, group, input, prefix = "b",
+                         suffix = "") {
+  if ((prefix == "i") && (suffix == "")) {
+    suffix <- c("1", "2")
+  }
+  if (is.null(group) || (group == "")) return()
+
+  for (suff in suffix) {
+    sc <- sprintf("%s_scenario%s", prefix, suff)
+    scenarios <- get_scenarios(touchstone, disease, group)
+    scenario <- update_dropdown_keep(session, sc, scenarios, input[[sc]])
+  }
+  update_scenario(session, touchstone, disease, group, scenario, input, prefix)
+}
+
+update_group_ts <- function(session, touchstone1, touchstone2, disease, group,
+                            input, prefix = "bts", suffix = "") {
+  if ((prefix == "its") && (suffix == "")) {
+    suffix <- c("1", "2")
+  }
+  if (is.null(group) || (group == "")) return()
+
+  for (suff in suffix) {
+    sc <- sprintf("%s_scenario%s", prefix, suff)
+    scenarios <- get_scenarios_ts(touchstone1, touchstone2, disease, group)
+    scenario <- update_dropdown_keep(session, sc, scenarios, input[[sc]])
+  }
+  update_scenario_ts(session, touchstone1, touchstone2, disease, group,
+                     scenario, input, prefix)
+}
+
+update_group_mg <- function(session, touchstone, disease, group1, group2,
+                            input, prefix = "bmg", suffix = "") {
+  if (is.null(group1) || is.null(group2) || (group1 == "") || (group2 == "")) return()
+  sc <- sprintf("%s_scenario", prefix)
+  scenarios <- get_scenarios_mg(touchstone, disease, group1, group2)
+  scenario <- update_dropdown_keep(session, sc, scenarios, input[[sc]])
+  update_scenario_mg(session, touchstone, disease, group1, group2,
+                     scenario, input, prefix)
+}
+
+update_scenario <- function(session, touchstone, disease, group, scenario, input,
+                            prefix = "b", extra_scenario = NULL) {
+  if (is.null(scenario) || (scenario == "")) return()
+  cc <- paste(prefix, "country", sep = "_")
+  countries <- get_countries(touchstone, disease, group, scenario, extra_scenario)
+  country <- update_dropdown_keep(session, cc, countries, input[[cc]])
+  update_country(session, touchstone, disease, group, scenario,
+                 country, input, prefix)
+}
+
+update_scenario_ts <- function(session, touchstone1, touchstone2, disease, group,
+                               scenario, input, prefix = "bts") {
+  if (is.null(scenario) || (scenario == "")) return()
+  cc <- paste(prefix, "country", sep = "_")
+  countries <- get_countries_ts(touchstone1, touchstone2, disease, group,
+                                scenario)
+  country <- update_dropdown_keep(session, cc, countries, input[[cc]])
+  update_country_ts(session, touchstone1, touchstone2, disease, group,
+                    scenario, country, input, prefix)
+}
+
+update_scenario_mg <- function(session, touchstone, disease, group1, group2,
+                               scenario, input, prefix = "bmg") {
+  if (is.null(scenario) || (scenario == "")) return()
+  if (is.null(group1) || is.null(group2) || (group1 == "") || (group2 == "")) return()
+  cc <- paste(prefix, "country", sep = "_")
+  countries <- get_countries_mg(touchstone, disease, group1, group2,
+                                scenario)
+  country <- update_dropdown_keep(session, cc, countries, input[[cc]])
+  update_country_mg(session, touchstone, disease, group1, group2,
+                    scenario, country, input, prefix)
+}
+
+update_country <- function(session, touchstone, disease, group, scenario, country,
+                           input, prefix = "b") {
+  if (is.null(country) || (country == "")) return()
+  o <- paste(prefix, "outcome", sep = "_")
+  outcomes <- get_outcomes(touchstone, disease, group, scenario, country)
+  update_dropdown_keep(session, o, outcomes, input[[o]])
+}
+
+update_country_ts <- function(session, touchstone1, touchstone2, disease, group,
+                              scenario, country, input, prefix = "bts") {
+  if (is.null(country) || (country == "")) return()
+  o <- paste(prefix, "outcome", sep = "_")
+  outcomes <- get_outcomes_ts(touchstone1, touchstone2, disease, group,
+                              scenario, country)
+  update_dropdown_keep(session, o, outcomes, input[[o]])
+}
+
+update_country_mg <- function(session, touchstone, disease, group1, group2,
+                              scenario, country, input, prefix = "bmg") {
+  if (is.null(country) || (country == "")) return()
+  if (is.null(group1) || is.null(group2) || (group1 == "") || (group2 == "")) return()
+  if (is.null(scenario) || (scenario == "")) return()
+  o <- paste(prefix, "outcome", sep = "_")
+  outcomes <- get_outcomes_mg(touchstone, disease, group1, group2,
+                              scenario, country)
+  update_dropdown_keep(session, o, outcomes, input[[o]])
+}

From 9062dd6cca27acec8f4f03daad766879debfe688 Mon Sep 17 00:00:00 2001
From: Wes Hinsley <w.hinsley@imperial.ac.uk>
Date: Thu, 26 Mar 2026 12:03:14 +0000
Subject: [PATCH 09/29] Add stoner::stochastic_explorer()

---
 NAMESPACE                  |  1 +
 R/stochastic_graphs.R      |  9 ++++++++-
 man/stochastic_explorer.Rd | 14 ++++++++++++++
 3 files changed, 23 insertions(+), 1 deletion(-)
 create mode 100644 man/stochastic_explorer.Rd

diff --git a/NAMESPACE b/NAMESPACE
index e07131a..e9fe275 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -1,6 +1,7 @@
 # Generated by roxygen2: do not edit by hand
 
 export(fetch_packit)
+export(stochastic_explorer)
 export(stone_dump)
 export(stone_extract)
 export(stone_load)
diff --git a/R/stochastic_graphs.R b/R/stochastic_graphs.R
index 1a6e2e2..b6723c4 100644
--- a/R/stochastic_graphs.R
+++ b/R/stochastic_graphs.R
@@ -153,6 +153,13 @@ prepare_central_data <- function(packit_id, packit_file,
   central
 }
 
-run_app <- function() {
+##' Launch a Shiny app to allow interactive plotting of
+##' standardised stochastic data, burden estimates,
+##' impacts, comparisons between touchstones, and
+##' comparisons between modelling groups.
+##'
+##' @export
+##' @title Stochastic plot
+stochastic_explorer <- function() {
   shiny::runApp(system.file("app", package = "stoner"))
 }
diff --git a/man/stochastic_explorer.Rd b/man/stochastic_explorer.Rd
new file mode 100644
index 0000000..5af6cd7
--- /dev/null
+++ b/man/stochastic_explorer.Rd
@@ -0,0 +1,14 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/stochastic_graphs.R
+\name{stochastic_explorer}
+\alias{stochastic_explorer}
+\title{Stochastic plot}
+\usage{
+stochastic_explorer()
+}
+\description{
+Launch a Shiny app to allow interactive plotting of
+standardised stochastic data, burden estimates,
+impacts, comparisons between touchstones, and
+comparisons between modelling groups.
+}

From 6e4121ccd8df619d3f338670946ccb398e9d0f84 Mon Sep 17 00:00:00 2001
From: Wes Hinsley <w.hinsley@imperial.ac.uk>
Date: Thu, 26 Mar 2026 20:34:44 +0000
Subject: [PATCH 10/29] Working app

---
 inst/app/app.R   | 621 ++++++++++++++++++++---------------------------
 inst/app/utils.R | 316 +++++++++---------------
 2 files changed, 377 insertions(+), 560 deletions(-)

diff --git a/inst/app/app.R b/inst/app/app.R
index 1bc368b..d1ed9e5 100644
--- a/inst/app/app.R
+++ b/inst/app/app.R
@@ -8,6 +8,63 @@ source("utils.R")
 
 app_ui <- function() {
   useShinyjs()
+
+  make_multi_input <- function(prefix, name, count, choices,
+                               n1 = name, n2 = NULL) {
+    if (count == 1) {
+      return(list(selectInput(sprintf("%s_%s", prefix, name),
+                              paste0(name, ":"),
+                              choices = choices, selected = NULL)))
+    } else {
+      return(list(selectInput(sprintf("%s_%s1", prefix, name), n1,
+                              choices = choices, selected = NULL),
+                  selectInput(sprintf("%s_%s2", prefix, name), n2,
+                              choices = choices, selected = NULL)))
+    }
+  }
+
+  make_panel <- function(title, prefix, n_touchstones, n_scenarios, n_groups) {
+    c0 <- character(0)
+    ts <- get_touchstones()
+
+    touchstones <- make_multi_input(prefix, "touchstone", n_touchstones, ts,
+                                    "Touchstone 1:", "Touchstone 2:")
+    scenarios <- make_multi_input(prefix, "scenario", n_scenarios, c0,
+                                    "Base Scenario:", "Compare Scenario:")
+    groups <- make_multi_input(prefix, "group", n_groups, c0,
+                               "Group 1:", "Group 2:")
+
+    if ((n_touchstones > 1) || (n_groups > 1)) {
+      graphs <- list(
+        plotOutput(sprintf("%s_main_plot1", prefix)),
+        plotOutput(sprintf("%s_main_plot2", prefix)))
+    } else {
+      graphs <- list(
+        plotOutput(sprintf("%s_main_plot", prefix)))
+    }
+
+    sidebar <- append(touchstones, list(
+      selectInput(sprintf("%s_disease", prefix), "Disease:", choices = c0)))
+    sidebar <- append(sidebar, groups)
+    sidebar <- append(sidebar, scenarios)
+    sidebar <- append(sidebar, list(
+      selectInput(sprintf("%s_country", prefix), "Country:", choices = c0),
+      selectInput(sprintf("%s_outcome", prefix), "Y-Axis:", choices = c0),
+      radioButtons(sprintf("%s_year", prefix), "X-Axis:", inline = TRUE,
+                   choices = c("Calendar", "Cohort")),
+      radioButtons(sprintf("%s_ages", prefix), "Age:", inline = TRUE,
+                   choices = c("All", "Under 5")),
+      checkboxGroupInput(
+        sprintf("%s_options", prefix), "Plot Options:",
+        choices = c("Quantiles", "Median", "Mean", "Log-Y"),
+        selected = c("Quantiles", "Median", "Mean", "Log-Y"),
+        inline = TRUE),
+      actionButton(sprintf("%s_plot_btn", prefix), "Plot")
+    ))
+
+    tabPanel(title, sidebarLayout(sidebarPanel(sidebar), mainPanel(graphs)))
+  }
+
   fluidPage(
     tags$style(HTML("
       .selectize-input {
@@ -51,377 +108,225 @@ app_ui <- function() {
     ")),
     titlePanel("Stochastic Explorer"),
     tabsetPanel(
-      tabPanel("Burden",
-        sidebarLayout(
-          sidebarPanel(
-            selectInput("b_touchstone", "Touchstone:",
-                        choices = get_touchstones(), selected = NULL),
-            selectInput("b_disease", "Disease:",
-                        choices = character(0)),
-            selectInput("b_group", "Group:",
-                        choices = character(0)),
-            selectInput("b_scenario", "Scenario:",
-                        choices = character(0)),
-            selectInput("b_country", "Country:",
-                        choices = character(0)),
-            selectInput("b_outcome", "Y-Axis:",
-                        choices = character(0)),
-            radioButtons("b_year", "X-Axis:", inline = TRUE,
-                         choices = c("Calendar", "Cohort")),
-            radioButtons("b_ages", "Age:", inline = TRUE,
-                         choices = c("All", "Under 5")),
-            checkboxGroupInput(
-              "b_options", "Plot Options:",
-              choices = c("Quantiles", "Median", "Mean", "Log-Y"),
-              selected = c("Quantiles", "Median", "Mean", "Log-Y"),
-              inline = TRUE),
-            actionButton("b_plot_btn", "Plot")
-          ), mainPanel(
-            plotOutput("b_main_plot")
-          )
-        )
-      ),
-      tabPanel("Burden/TS",
-        sidebarLayout(
-          sidebarPanel(
-            selectInput("bts_touchstone1", "Touchstone 1:",
-              choices = get_touchstones(), selected = NULL),
-            selectInput("bts_touchstone2", "Touchstone 2:",
-              choices = get_touchstones(), selected = NULL),
-            selectInput("bts_disease", "Disease:",
-              choices = character(0)),
-            selectInput("bts_group", "Group:",
-              choices = character(0)),
-            selectInput("bts_scenario", "Scenario:",
-              choices = character(0)),
-            selectInput("bts_country", "Country:",
-              choices = character(0)),
-            selectInput("bts_outcome", "Y-Axis:",
-              choices = character(0)),
-            radioButtons("bts_year", "X-Axis:", inline = TRUE,
-              choices = c("Calendar", "Cohort")),
-            radioButtons("bts_ages", "Age:", inline = TRUE,
-              choices = c("All", "Under 5")),
-            checkboxGroupInput("bts_options", "Plot Options:",
-              choices = c("Quantiles", "Median", "Mean", "Log-Y"),
-                       selected = c("Quantiles", "Median", "Mean", "Log-Y"),
-                       inline = TRUE),
-            actionButton("bts_plot_btn", "Plot")
-          ), mainPanel(
-            plotOutput("bts_main_plot1"),
-            plotOutput("bts_main_plot2")
-          )
-        )
-      ),
-      tabPanel("Burden/MG",
-        sidebarLayout(
-          sidebarPanel(
-            selectInput("bmg_touchstone", "Touchstone:",
-                        choices = get_touchstones(), selected = NULL),
-            selectInput("bmg_disease", "Disease:",
-                        choices = character(0)),
-            selectInput("bmg_group1", "Group 1:",
-                        choices = character(0)),
-            selectInput("bmg_group2", "Group 2:",
-                        choices = character(0)),
-            selectInput("bmg_scenario", "Scenario:",
-                        choices = character(0)),
-            selectInput("bmg_country", "Country:",
-                        choices = character(0)),
-            selectInput("bmg_outcome", "Y-Axis:",
-                        choices = character(0)),
-            radioButtons("bmg_year", "X-Axis:", inline = TRUE,
-                         choices = c("Calendar", "Cohort")),
-            radioButtons("bmg_ages", "Age:", inline = TRUE,
-                         choices = c("All", "Under 5")),
-            checkboxGroupInput("bmg_options", "Plot Options:",
-                         choices = c("Quantiles", "Median", "Mean", "Log-Y"),
-                         selected = c("Quantiles", "Median", "Mean", "Log-Y"),
-                         inline = TRUE),
-            actionButton("bmg_plot_btn", "Plot")
-          ), mainPanel(
-            plotOutput("bmg_main_plot1"),
-            plotOutput("bmg_main_plot2")
-          )
-        )
-      ),
-      tabPanel("Impact",
-        sidebarLayout(
-          sidebarPanel(
-            selectInput("i_touchstone", "Touchstone:",
-                       choices = get_touchstones(), selected = NULL),
-            selectInput("i_disease", "Disease:",
-                       choices = character(0)),
-            selectInput("i_group", "Group:",
-                       choices = character(0)),
-            selectInput("i_scenario1", "Base Scenario:",
-                       choices = character(0)),
-            selectInput("i_scenario2", "Compare Scenario:",
-                       choices = character(0)),
-            selectInput("i_country", "Country:",
-                       choices = character(0)),
-            selectInput("i_outcome", "Y-Axis:",
-                       choices = character(0)),
-            radioButtons("i_year", "X-Axis:", inline = TRUE,
-                        choices = c("Calendar", "Cohort")),
-            radioButtons("i_ages", "Age:", inline = TRUE,
-                        choices = c("All", "Under 5")),
-            checkboxGroupInput(
-                 "i_options", "Plot Options:",
-                 choices = c("Quantiles", "Median", "Mean", "Log-Y"),
-                 selected = c("Quantiles", "Median", "Mean", "Log-Y"),
-                 inline = TRUE),
-            actionButton("i_plot_btn", "Plot")
-          ), mainPanel(
-           plotOutput("i_main_plot")
-          )
-        )
-      )
+      make_panel("Burden", "b", 1, 1, 1),
+      make_panel("Burden/TS", "bts", 2, 1, 1),
+      make_panel("Burden/MG", "bmg", 1, 1, 2),
+      make_panel("Impact", "i", 1, 2, 1),
+      make_panel("Impact/TS", "its", 2, 2, 1),
+      make_panel("Impact/MG", "img", 1, 2, 2)
     )
   )
 }
 
 app_server <- function(input, output, session) {
 
-  # Events for the basic burden tab
-
-  observeEvent(input$b_touchstone, {
-    update_touchstone(session, input$b_touchstone, input) })
-  observeEvent(input$b_disease, {
-    update_disease(session, input$b_touchstone, input$b_disease, input) })
-  observeEvent(input$b_group, {
-    update_group(session, input$b_touchstone, input$b_disease, input$b_group, input) })
-  observeEvent(input$b_scenario, {
-    update_scenario(session, input$b_touchstone, input$b_disease, input$b_group,
-                    input$b_scenario, input) })
-  observeEvent(input$b_country, {
-    update_country(session, input$b_touchstone, input$b_disease, input$b_group,
-                   input$b_scenario, input$b_country, input) })
-
-
-  plot_obj_b <- eventReactive(input$b_plot_btn, {
-    ages <- NULL
-    if (input$b_ages == "Under 5") {
-      ages <- 0:4
-    }
-    stone_stochastic_graph(
-      base = data_dir,
-      touchstone = input$b_touchstone,
-      disease = input$b_disease,
-      group = input$b_group,
-      country = input$b_country,
-      scenario = input$b_scenario,
-      outcome = input$b_outcome,
-      by_cohort = input$b_year == "Cohort",
-      ages = ages,
-      log = "Log-Y" %in% input$b_options,
-      include_median = "Median" %in% input$b_options,
-      include_quantiles = "Quantiles" %in% input$b_options,
-      include_mean = "Mean" %in% input$b_options,
-    )
-  })
-
-  output$b_main_plot <- renderPlot({
-    plot_obj_b()
-  })
-
-  # Events for the impact tab
-
-  observeEvent(input$i_touchstone, {
-    update_touchstone(session, input$i_touchstone, input, "i") })
-  observeEvent(input$i_disease, {
-    update_disease(session, input$i_touchstone, input$i_disease, input, "i") })
-  observeEvent(input$i_group, {
-    update_group(session, input$i_touchstone, input$i_disease,
-                 input$i_group, input, "i", "1")
-    update_group(session, input$i_touchstone, input$i_disease,
-                 input$i_group, input, "i", "2")
-  })
-  observeEvent(input$i_scenario1, {
-    update_scenario(session, input$i_touchstone, input$i_disease, input$i_group,
-                    input$i_scenario1, input, "i",
-                    extra_scenario = input$i_scenario2) })
-  observeEvent(input$i_scenario2, {
-    update_scenario(session, input$i_touchstone, input$i_disease, input$i_group,
-                    input$i_scenario2, input, "i",
-                    extra_scenario = input$i_scenario2) })
-  observeEvent(input$i_country, {
-    update_country(session, input$i_touchstone, input$i_disease, input$i_group,
-                   input$i_scenario1, input$i_country, input, "i") })
-
-  plot_obj_i <- eventReactive(input$i_plot_btn, {
-
-    if (input$i_scenario1 == input$i_scenario2) {
-      showModal(modalDialog(
-        title = "Warning",
-        "Same scenario selected - impact will be zero.",
-        easyClose = TRUE,
-        footer = modalButton("OK")
-      ))
-      return(NULL)
+  add_observers <- function(prefix, n_touchstone = 1, n_group = 1,
+                            n_scenario = 1) {
+    if (n_touchstone == 1) {
+      it1 <- sprintf("%s_touchstone", prefix)
+      it2 <- NULL
+      observeEvent(input[[it1]], {
+                   update_touchstone(session,
+                                     input[[it1]],
+                                     NULL,
+                                     input,
+                                     prefix)})
+
+    } else {
+      it1 <- sprintf("%s_touchstone1", prefix)
+      it2 <- sprintf("%s_touchstone2", prefix)
+
+      observeEvent(input[[it1]], {
+                   update_touchstone(session,
+                                      input[[it1]],
+                                      input[[it2]],
+                                      input,
+                                      prefix)})
+
+      observeEvent(input[[it2]], {
+                   update_touchstone(session,
+                                     input[[it1]],
+                                     input[[it2]],
+                                     input,
+                                     prefix)})
     }
 
-    ages <- NULL
-    if (input$i_ages == "Under 5") {
-      ages <- 0:4
-    }
-    stone_stochastic_graph(
-      base = data_dir,
-      touchstone = input$i_touchstone,
-      disease = input$i_disease,
-      group = input$i_group,
-      country = input$i_country,
-      scenario = input$i_scenario1,
-      scenario2 = input$i_scenario2,
-      outcome = input$i_outcome,
-      by_cohort = input$i_year == "Cohort",
-      ages = ages,
-      log = "Log-Y" %in% input$i_options,
-      include_median = "Median" %in% input$i_options,
-      include_quantiles = "Quantiles" %in% input$i_options,
-      include_mean = "Mean" %in% input$i_options,
-    )
-  })
-
-  output$i_main_plot <- renderPlot({
-    plot_obj_i()
-  })
-
-  # Events for comparing burdens on different touchstones
-
-  observeEvent(input$bts_touchstone1, {
-    update_touchstone_ts(session, input$bts_touchstone1, input$bts_touchstone2,
-                         input, "bts") })
-  observeEvent(input$bts_touchstone2, {
-    update_touchstone_ts(session, input$bts_touchstone1, input$bts_touchstone2,
-                         input, "bts") })
-  observeEvent(input$bts_disease, {
-    update_disease_ts(session, input$bts_touchstone1, input$bts_touchstone2,
-                      input$bts_disease, input, "bts") })
-  observeEvent(input$bts_group, {
-    update_group_ts(session, input$bts_touchstone1, input$bts_touchstone2,
-                    input$bts_disease, input$bts_group, input, "bts")
-  })
-  observeEvent(input$bts_scenario, {
-    update_scenario_ts(session, input$bts_touchstone1, input$bts_touchstone2,
-                       input$bts_disease, input$bts_group, input$bts_scenario,
-                       input, "bts") })
-  observeEvent(input$bts_country, {
-    update_country_ts(session, input$bts_touchstone1, input$bts_touchstone2,
-                      input$bts_disease, input$bts_group, input$bts_scenario,
-                      input$bts_country, input, "bts") })
-
-  plot_obj_bts <- eventReactive(input$bts_plot_btn, {
-    if (input$bts_touchstone1 == input$bts_touchstone2) {
-      showModal(modalDialog(
-        title = "Warning",
-        "Same touchstones selected.",
-        easyClose = TRUE,
-        footer = modalButton("OK")
-      ))
-      return(NULL)
-    }
-    ages <- NULL
-    if (input$i_ages == "Under 5") {
-      ages <- 0:4
-    }
-    plots <- list()
-    for (ts in 1:2) {
-      plots[[ts]] <-
-        stone_stochastic_graph(
-          base = data_dir,
-          touchstone = if (ts == 1) input$bts_touchstone1 else input$bts_touchstone2,
-          disease = input$bts_disease,
-          group = input$bts_group,
-          country = input$bts_country,
-          scenario = input$bts_scenario,
-          outcome = input$bts_outcome,
-          by_cohort = input$bts_year == "Cohort",
-          ages = ages,
-          log = "Log-Y" %in% input$bts_options,
-          include_median = "Median" %in% input$bts_options,
-          include_quantiles = "Quantiles" %in% input$bts_options,
-          include_mean = "Mean" %in% input$bts_options,
-        )
-    }
-    list(plot1 = plots[[1]], plot2 = plots[[2]])
-  })
-
-  output$bts_main_plot1 <- renderPlot({
-    replayPlot(plot_obj_bts()$plot1)
-  })
-
-  output$bts_main_plot2 <- renderPlot({
-    replayPlot(plot_obj_bts()$plot2)
-  })
-
-  # Compare burden estimates for different groups
-
-  observeEvent(input$bmg_touchstone, {
-    update_touchstone(session, input$bmg_touchstone, input, "bmg") })
-  observeEvent(input$bmg_disease, {
-    update_disease_mg(session, input$bmg_touchstone, input$bmg_disease,
-                   input, "bmg") })
-  observeEvent(input$bmg_group1, {
-    update_group_mg(session, input$bmg_touchstone, input$bmg_disease,
-                    input$bmg_group1, input$bmg_group2, input, "bmg")
-  })
-  observeEvent(input$bmg_group2, {
-    update_group_mg(session, input$bmg_touchstone, input$bmg_disease,
-                    input$bmg_group1, input$bmg_group2, input, "bmg")
-  })
-  observeEvent(input$bmg_scenario, {
-    update_scenario_mg(session, input$bmg_touchstone, input$bmg_disease,
-                       input$bmg_group1, input$bmg_group2, input$bmg_scenario,
-                       input, "bmg") })
-  observeEvent(input$bmg_country, {
-    update_country_mg(session, input$bmg_touchstone, input$bmg_disease,
-                      input$bmg_group1, input$bmg_group2, input$bmg_scenario,
-                      input$bmg_country, input, "bmg") })
-
-  plot_obj_bmg <- eventReactive(input$bmg_plot_btn, {
-    if (input$bmg_group1 == input$bmg_group2) {
-      showModal(modalDialog(
-        title = "Warning",
-        "Same groups selected.",
-        easyClose = TRUE,
-        footer = modalButton("OK")
-      ))
-      return(NULL)
+    id <- sprintf("%s_disease", prefix)
+    observeEvent(input[[id]], {
+                 update_disease(session,
+                                input[[it1]],
+                                if (is.null(it2)) NULL else input[[it2]],
+                                input[[id]],
+                                input,
+                                prefix)})
+
+    if (n_group == 1) {
+      ig1 <- sprintf("%s_group", prefix)
+      ig2 <- NULL
+      observeEvent(input[[ig1]], {
+                   update_group(session,
+                                input[[it1]],
+                                if (is.null(it2)) NULL else input[[it2]],
+                                input[[id]],
+                                input[[ig1]],
+                                NULL,
+                                input,
+                                prefix)})
+
+    } else {
+      ig1 <- sprintf("%s_group1", prefix)
+      ig2 <- sprintf("%s_group2", prefix)
+      observeEvent(input[[ig1]], {
+                   update_group(session,
+                                input[[it1]],
+                                if (is.null(it2)) NULL else input[[it2]],
+                                input[[id]],
+                                input[[ig1]],
+                                input[[ig2]],
+                                input,
+                                prefix)})
+
+      observeEvent(input[[ig2]], {
+                   update_group(session,
+                                input[[it1]],
+                                if (is.null(it2)) NULL else input[[it2]],
+                                input[[id]],
+                                input[[ig1]],
+                                input[[ig2]],
+                                input,
+                                prefix)})
     }
-    ages <- NULL
-    if (input$i_ages == "Under 5") {
-      ages <- 0:4
+
+    if (n_scenario == 1) {
+      is1 <- sprintf("%s_scenario", prefix)
+      is2 <- NULL
+      observeEvent(input[[is1]], {
+                   update_scenario(session,
+                                   input[[it1]],
+                                   if (is.null(it2)) NULL else input[[it2]],
+                                   input[[id]],
+                                   input[[ig1]],
+                                   if (is.null(ig2)) NULL else input[[ig2]],
+                                   input[[is1]],
+                                   NULL,
+                                   input, prefix)})
+    } else {
+      is1 <- sprintf("%s_scenario1", prefix)
+      is2 <- sprintf("%s_scenario2", prefix)
+      observeEvent(input[[is1]], {
+                   update_scenario(session,
+                                   input[[it1]],
+                                   if (is.null(it2)) NULL else input[[it2]],
+                                   input[[id]],
+                                   input[[ig1]],
+                                   if (is.null(ig2)) NULL else input[[ig2]],
+                                   input[[is1]],
+                                   input[[is2]],
+                                   input,
+                                   prefix)})
+
+      observeEvent(input[[is2]], {
+                   update_scenario(session,
+                                   input[[it1]],
+                                   if (is.null(it2)) NULL else input[[it2]],
+                                   input[[id]],
+                                   input[[ig1]],
+                                   if (is.null(ig2)) NULL else input[[ig2]],
+                                   input[[is1]],
+                                   input[[is2]],
+                                   input,
+                                   prefix)})
     }
-    plots <- list()
-    for (ts in 1:2) {
-      plots[[ts]] <-
-        stone_stochastic_graph(
-          base = data_dir,
-          touchstone = input$bmg_touchstone,
-          disease = input$bmg_disease,
-          group = if (ts == 1) input$bmg_group1 else input$bmg_group2,
-          country = input$bmg_country,
-          scenario = input$bmg_scenario,
-          outcome = input$bmg_outcome,
-          by_cohort = input$bmg_year == "Cohort",
-          ages = ages,
-          log = "Log-Y" %in% input$bmg_options,
-          include_median = "Median" %in% input$bmg_options,
-          include_quantiles = "Quantiles" %in% input$bmg_options,
-          include_mean = "Mean" %in% input$bmg_options,
-        )
+
+    ic <- sprintf("%s_country", prefix)
+    observeEvent(input[[ic]], {
+                 update_country(session,
+                                input[[it1]],
+                                if (is.null(it2)) NULL else input[[it2]],
+                                input[[id]],
+                                input[[ig1]],
+                                if (is.null(ig2)) NULL else input[[ig2]],
+                                input[[is1]],
+                                if (is.null(is2)) NULL else input[[is2]],
+                                input[[ic]],
+                                input,
+                                prefix)})
+
+    n_graphs <- 1 + ((n_touchstone * n_group) > 1)
+    if (n_graphs == 1) {
+      graphs <- sprintf("%s_main_plot", prefix)
+    } else {
+      graphs <- sprintf("%s_main_plot%d", prefix, 1:2)
     }
-    list(plot1 = plots[[1]], plot2 = plots[[2]])
-  })
+    button <- sprintf("%s_plot_btn", prefix)
 
-  output$bmg_main_plot1 <- renderPlot({
-    replayPlot(plot_obj_bmg()$plot1)
-  })
+    plot_reactive <- eventReactive(input[[button]], {
+      ages <- NULL
+      if (input[[sprintf("%s_ages", prefix)]] == "Under 5") {
+        ages <- 0:4
+      }
 
-  output$bmg_main_plot2 <- renderPlot({
-    replayPlot(plot_obj_bmg()$plot2)
-  })
+      check <- function(pre, x1, x2, type) {
+        if (is.null(x2)) return(TRUE)
+        if ((grepl(pre, prefix)) && (input[[x1]] == input[[x2]])) {
+          showModal(modalDialog(
+            title = "Warning",
+            sprintf("Same %s selected. Dull comparison.", type),
+            easyClose = TRUE,
+            footer = modalButton("OK")
+          ))
+          return(FALSE)
+        }
+        TRUE
+      }
+
+      if (!check("mg", ig1, ig2, "group")) return(NULL)
+      if (!check("ts", it1, it2, "touchstone")) return(NULL)
+      if (!check("i", is1, is2, "scenario")) return(NULL)
+
+      list(
+        ages = ages,
+        opts = input[[sprintf("%s_options", prefix)]]
+      )
+    })
+
+    for (g in seq_along(graphs)) {
+      local({
+        gg <- g
+
+        output[[graphs[gg]]] <- renderPlot({
+          pr <- plot_reactive()
+          req(pr)
+
+          it <- if (n_touchstone == 1) it1 else if (gg == 1) it1 else it2
+          ig <- if (n_group == 1) ig1 else if (gg == 1) ig1 else ig2
+
+          stone_stochastic_graph(
+            base = data_dir,
+            touchstone = input[[it]],
+            disease = input[[id]],
+            group = input[[ig]],
+            country = input[[ic]],
+            scenario = input[[is1]],
+            scenario2 = if (is.null(is2)) NULL else input[[is2]],
+            outcome = input[[sprintf("%s_outcome", prefix)]],
+            by_cohort = input[[sprintf("%s_year", prefix)]] == "Cohort",
+            ages = pr$ages,
+            log = "Log-Y" %in% pr$opts,
+            include_median = "Median" %in% pr$opts,
+            include_quantiles = "Quantiles" %in% pr$opts,
+            include_mean = "Mean" %in% pr$opts
+          )
+        })
+      })
+    }
+  }
+  add_observers("b")
+  add_observers("bts", 2, 1, 1)
+  add_observers("bmg", 1, 2, 1)
+  add_observers("i", 1, 1, 2)
+  add_observers("its", 2, 1, 2)
+  add_observers("img", 1, 2, 2)
 }
 
 shinyApp(app_ui, app_server)
diff --git a/inst/app/utils.R b/inst/app/utils.R
index f68393f..1500413 100644
--- a/inst/app/utils.R
+++ b/inst/app/utils.R
@@ -9,101 +9,98 @@
 
 data_dir <- "//wpia-hn2.hpc.dide.ic.ac.uk/vimc_stochastics"
 
-
 get_touchstones <- function() {
-  basename(list.dirs(data_dir, recursive = FALSE))
-}
-
-get_diseases <- function(touchstone) {
-  path <- file.path(data_dir, touchstone)
-  dirs <- basename(list.dirs(path, recursive = FALSE))
-  unique(unlist(lapply(strsplit(dirs, "_"), `[[`, 1)))
-}
-
-get_diseases_ts <- function(touchstone1, touchstone2) {
-  res <- get_diseases(touchstone1)
-  res[res %in% get_diseases(touchstone2)]
+  sort(unique(basename(list.dirs(data_dir, recursive = FALSE))))
 }
 
-get_groups <- function(touchstone, disease) {
-  path <- file.path(data_dir, touchstone)
-  dirs <- basename(list.dirs(path, recursive = FALSE))
-  dirs <- dirs[substr(dirs, 1, nchar(disease) + 1) == paste0(disease, "_")]
-  unique(substring(dirs, nchar(disease) + 2))
+get_diseases <- function(touchstone1, touchstone2) {
+  lookup <- function(touchstone, res = NULL) {
+    if (is.null(touchstone)) return(res)
+    path <- file.path(data_dir, touchstone)
+    dirs <- basename(list.dirs(path, recursive = FALSE))
+    unique(unlist(lapply(strsplit(dirs, "_"), `[[`, 1)))
+  }
+  res <- lookup(touchstone1)
+  sort(unique(res[res %in% lookup(touchstone2, res)]))
 }
 
-get_groups_ts <- function(touchstone1, touchstone2, disease) {
-  res <- get_groups(touchstone1, disease)
-  res[res %in% get_groups(touchstone2, disease)]
+get_groups <- function(touchstone1, touchstone2, disease) {
+  lookup <- function(touchstone, disease, res = NULL) {
+    if (is.null(touchstone)) return(res)
+    path <- file.path(data_dir, touchstone)
+    dirs <- basename(list.dirs(path, recursive = FALSE))
+    dirs <- dirs[substr(dirs, 1, nchar(disease) + 1) == paste0(disease, "_")]
+    unique(substring(dirs, nchar(disease) + 2))
+  }
+  res <- lookup(touchstone1, disease)
+  sort(unique(res[res %in% lookup(touchstone2, disease, res)]))
 }
 
-get_scenarios <- function(touchstone, disease, group) {
-  path <- file.path(data_dir, touchstone, paste(disease, group, sep = "_"))
-  files <- basename(list.files(path, recursive = FALSE))
-  files <- gsub(".pq", "", files)
-  files <- substring(files, nchar(group) + 2)
-  ends <- unlist(lapply(gregexpr("_", files), `[[`, 1)) - 1
-  unique(substring(files, 1, ends))
-}
+get_scenarios <- function(touchstone1, touchstone2, disease, group1, group2) {
+  lookup <- function(touchstone, disease, group, res = NULL) {
+    if ((is.null(touchstone)) || (is.null(group))) return(res)
+    path <- file.path(data_dir, touchstone, paste(disease, group, sep = "_"))
+    files <- basename(list.files(path, recursive = FALSE))
+    files <- gsub(".pq", "", files)
+    files <- substring(files, nchar(group) + 2)
+    ends <- unlist(lapply(gregexpr("_", files), `[[`, 1)) - 1
+    unique(substring(files, 1, ends))
+  }
 
-get_scenarios_ts <- function(touchstone1, touchstone2, disease, group) {
-  res <- get_scenarios(touchstone1, disease, group)
-  res[res %in% get_scenarios(touchstone2, disease, group)]
+  res <- lookup(touchstone1, disease, group1)
+  res <- res[res %in% lookup(touchstone2, disease, group1, res)]
+  res <- res[res %in% lookup(touchstone1, disease, group2, res)]
+  sort(unique(res[res %in% lookup(touchstone2, disease, group2, res)]))
 }
 
-get_scenarios_mg <- function(touchstone, disease, group1, group2) {
-  res <- get_scenarios(touchstone, disease, group1)
-  res[res %in% get_scenarios(touchstone, disease, group2)]
-}
+get_countries <- function(touchstone1, touchstone2, disease, group1, group2,
+                          scenario1, scenario2) {
+  lookup <- function(touchstone, disease, group, scenario, res = NULL) {
+    if ((is.null(touchstone)) || (is.null(group)) || (is.null(scenario))) {
+      return(res)
+    }
 
-get_countries <- function(touchstone, disease, group, scenario, extra_scenario = NULL) {
-  scenarios <- c(scenario, extra_scenario)
-  final <- NULL
-  for (scenario in scenarios) {
     path <- file.path(data_dir, touchstone, paste(disease, group, sep = "_"))
     files <- basename(list.files(path, recursive = FALSE))
     files <- gsub(".pq", "", files)
     files <- substring(files, nchar(group) + 2)
     files <- files[substring(files, 1, nchar(scenario) + 1) == paste0(scenario, "_")]
     ends <- unlist(lapply(gregexpr("_", files), `[[`, 1)) + 1
-    if (is.null(final)) final <- unique(substring(files, ends))
-    else final <- final[final %in% unique(substring(files, ends))]
+    unique(substring(files, ends))
   }
-  final
-}
-
-get_countries_ts <- function(touchstone1, touchstone2, disease, group, scenario) {
-  res <- get_countries(touchstone1, disease, group, scenario)
-  res[res %in% get_countries(touchstone2, disease, group, scenario)]
-}
-
-get_countries_mg <- function(touchstone, disease, group1, group2, scenario) {
-  res <- get_countries(touchstone, disease, group1, scenario)
-  res[res %in% get_countries(touchstone, disease, group2, scenario)]
-}
-
-get_outcomes <- function(touchstone, disease, group, scenario, country) {
-  path <- file.path(data_dir, touchstone, paste(disease, group, sep = "_"))
-  thefile <- sprintf("%s_%s_%s.pq", group, scenario, country)
-  tbl <- arrow::read_parquet(file.path(path, thefile), col_select = NULL)
-  cols <- names(tbl)
-  sort(cols[!cols %in% c("disease", "run_id", "year", "age",
-                    "country", "cohort_size")])
-}
-
-get_outcomes_ts <- function(touchstone1, touchstone2, disease, group,
-                            scenario, country) {
-  res <- get_outcomes(touchstone1, disease, group, scenario, country)
-  res[res %in% get_outcomes(touchstone2, disease, group, scenario, country)]
-}
-
-get_outcomes_mg <- function(touchstone, disease, group1, group2, scenario,
-                            country) {
-  res <- get_outcomes(touchstone, disease, group1, scenario, country)
-  res[res %in% get_outcomes(touchstone, disease, group2, scenario, country)]
+  res <- lookup(touchstone1, disease, group1, scenario1)
+  res <- res[res %in% lookup(touchstone2, disease, group1, scenario1, res)]
+  res <- res[res %in% lookup(touchstone1, disease, group2, scenario1, res)]
+  res <- res[res %in% lookup(touchstone2, disease, group2, scenario1, res)]
+  res <- res[res %in% lookup(touchstone1, disease, group1, scenario2, res)]
+  res <- res[res %in% lookup(touchstone2, disease, group1, scenario2, res)]
+  res <- res[res %in% lookup(touchstone1, disease, group2, scenario2, res)]
+  sort(unique(res[res %in% lookup(touchstone2, disease, group2, scenario2, res)]))
+}
+
+get_outcomes <- function(touchstone1, touchstone2, disease, group1, group2,
+                         scenario1, scenario2, country) {
+  lookup <- function(touchstone, disease, group, scenario, country, res = NULL) {
+    if ((is.null(touchstone)) || (is.null(group)) || (is.null(scenario))) {
+      return(res)
+    }
+    path <- file.path(data_dir, touchstone, paste(disease, group, sep = "_"))
+    thefile <- sprintf("%s_%s_%s.pq", group, scenario, country)
+    tbl <- arrow::read_parquet(file.path(path, thefile), col_select = NULL)
+    cols <- names(tbl)
+    sort(cols[!cols %in% c("disease", "run_id", "year", "age",
+                           "country", "cohort_size")])
+  }
+  res <- lookup(touchstone1, disease, group1, scenario1, country)
+  res <- res[res %in% lookup(touchstone2, disease, group1, scenario1, country, res)]
+  res <- res[res %in% lookup(touchstone1, disease, group2, scenario1, country, res)]
+  res <- res[res %in% lookup(touchstone2, disease, group2, scenario1, country, res)]
+  res <- res[res %in% lookup(touchstone1, disease, group1, scenario2, country, res)]
+  res <- res[res %in% lookup(touchstone2, disease, group1, scenario2, country, res)]
+  res <- res[res %in% lookup(touchstone1, disease, group2, scenario2, country, res)]
+  sort(unique(res[res %in% lookup(touchstone2, disease, group2, scenario2, country, res)]))
 }
 
-
 # GUI helpers.
 
 # Update a dropdown. Keep the previously selected value if possible,
@@ -122,149 +119,64 @@ update_dropdown_keep <- function(session, id, choices, selected) {
 # "bts" for multi-touchstone burden, "its" for multi-touchstone impact
 # "bmg" for multi-group burden, and "img" for multi-group impact.
 
-update_touchstone <- function(session, touchstone, input, prefix = "b") {
-  d <- paste(prefix, "disease", sep = "_")
-  diseases <- get_diseases(touchstone)
-  disease <- update_dropdown_keep(session, d, diseases, input[[d]])
-  update_disease(session, touchstone, disease, input, prefix)
-}
+# It's messy as sometimes we have two touchstones, groups, or scenarios.
 
-update_touchstone_ts <- function(session, touchstone1, touchstone2,
-                                 input, prefix = "bts") {
+update_touchstone <- function(session, ts1, ts2, input, prefix) {
   d <- paste(prefix, "disease", sep = "_")
-  diseases <- get_diseases_ts(touchstone1, touchstone2)
+  diseases <- get_diseases(ts1, ts2)
   disease <- update_dropdown_keep(session, d, diseases, input[[d]])
-  update_disease_ts(session, touchstone1, touchstone2, disease, input, prefix)
+  update_disease(session, ts1, ts2, disease, input, prefix)
 }
 
-update_disease <- function(session, touchstone, disease, input, prefix = "b") {
+update_disease <- function(session, ts1, ts2, disease, input, prefix) {
   if (is.null(disease) || (disease == "")) return()
-  g <- paste(prefix, "group", sep = "_")
-  groups <- get_groups(touchstone, disease)
-  group <- update_dropdown_keep(session, g, groups, input[[g]])
-  update_group(session, touchstone, disease, group, input, prefix)
-}
-
-update_disease_ts <- function(session, touchstone1, touchstone2, disease,
-                              input, prefix = "bts") {
-  if (is.null(disease) || (disease == "")) return()
-  g <- paste(prefix, "group", sep = "_")
-  groups <- get_groups_ts(touchstone1, touchstone2, disease)
-  group <- update_dropdown_keep(session, g, groups, input[[g]])
-  update_group_ts(session, touchstone1, touchstone2, disease, group,
-                  input, prefix)
-}
-
-update_disease_mg <- function(session, touchstone, disease, input,
-                              prefix = "bmg") {
-  if (is.null(disease) || (disease == "")) return()
-  g1 <- paste(prefix, "group1", sep = "_")
-  g2 <- paste(prefix, "group2", sep = "_")
-  groups <- get_groups(touchstone, disease)
-  group1 <- update_dropdown_keep(session, g1, groups, input[[g1]])
-  group2 <- update_dropdown_keep(session, g2, groups, input[[g2]])
-  update_group_mg(session, touchstone, disease, group1, group2,
-                  input, prefix)
-}
-
-update_group <- function(session, touchstone, disease, group, input, prefix = "b",
-                         suffix = "") {
-  if ((prefix == "i") && (suffix == "")) {
-    suffix <- c("1", "2")
-  }
-  if (is.null(group) || (group == "")) return()
-
-  for (suff in suffix) {
-    sc <- sprintf("%s_scenario%s", prefix, suff)
-    scenarios <- get_scenarios(touchstone, disease, group)
-    scenario <- update_dropdown_keep(session, sc, scenarios, input[[sc]])
+  groups <- get_groups(ts1, ts2, disease)
+  if (grepl("mg", prefix)) {
+    g1 <- paste(prefix, "group1", sep = "_")
+    g2 <- paste(prefix, "group2", sep = "_")
+    group1 <- update_dropdown_keep(session, g1, groups, input[[g1]])
+    group2 <- update_dropdown_keep(session, g2, groups, input[[g2]])
+    update_group(session, ts1, ts2, disease, group1, group2, input, prefix)
+  } else {
+    g <- paste(prefix, "group", sep = "_")
+    group1 <- update_dropdown_keep(session, g, groups, input[[g]])
+    group2 <- NULL
   }
-  update_scenario(session, touchstone, disease, group, scenario, input, prefix)
-}
-
-update_group_ts <- function(session, touchstone1, touchstone2, disease, group,
-                            input, prefix = "bts", suffix = "") {
-  if ((prefix == "its") && (suffix == "")) {
-    suffix <- c("1", "2")
-  }
-  if (is.null(group) || (group == "")) return()
-
-  for (suff in suffix) {
-    sc <- sprintf("%s_scenario%s", prefix, suff)
-    scenarios <- get_scenarios_ts(touchstone1, touchstone2, disease, group)
-    scenario <- update_dropdown_keep(session, sc, scenarios, input[[sc]])
+  update_group(session, ts1, ts2, disease, group1, group2, input, prefix)
+}
+
+update_group <- function(session, ts1, ts2, disease, g1, g2, input, prefix) {
+  if (is.null(g1) || (g1 == "")) return()
+  scenarios <- get_scenarios(ts1, ts2, disease, g1, g2)
+  if (grepl("i", prefix)) {
+    sc <- sprintf("%s_scenario1", prefix)
+    s1 <- update_dropdown_keep(session, sc, scenarios, input[[sc]])
+    sc <- sprintf("%s_scenario2", prefix)
+    s2 <- update_dropdown_keep(session, sc, scenarios, input[[sc]])
+  } else {
+    sc <- sprintf("%s_scenario", prefix)
+    s1 <- update_dropdown_keep(session, sc, scenarios, input[[sc]])
+    s2 <- NULL
   }
-  update_scenario_ts(session, touchstone1, touchstone2, disease, group,
-                     scenario, input, prefix)
+  update_scenario(session, ts1, ts2, disease, g1, g2, s1, s2, input, prefix)
 }
 
-update_group_mg <- function(session, touchstone, disease, group1, group2,
-                            input, prefix = "bmg", suffix = "") {
-  if (is.null(group1) || is.null(group2) || (group1 == "") || (group2 == "")) return()
-  sc <- sprintf("%s_scenario", prefix)
-  scenarios <- get_scenarios_mg(touchstone, disease, group1, group2)
-  scenario <- update_dropdown_keep(session, sc, scenarios, input[[sc]])
-  update_scenario_mg(session, touchstone, disease, group1, group2,
-                     scenario, input, prefix)
-}
+update_scenario <- function(session, ts1, ts2, disease, g1, g2, s1, s2,
+                            input, prefix) {
+  if (is.null(s1) || (s1 == "")) return()
 
-update_scenario <- function(session, touchstone, disease, group, scenario, input,
-                            prefix = "b", extra_scenario = NULL) {
-  if (is.null(scenario) || (scenario == "")) return()
   cc <- paste(prefix, "country", sep = "_")
-  countries <- get_countries(touchstone, disease, group, scenario, extra_scenario)
+  countries <- get_countries(ts1, ts2, disease, g1, g2, s1, s2)
   country <- update_dropdown_keep(session, cc, countries, input[[cc]])
-  update_country(session, touchstone, disease, group, scenario,
-                 country, input, prefix)
-}
-
-update_scenario_ts <- function(session, touchstone1, touchstone2, disease, group,
-                               scenario, input, prefix = "bts") {
-  if (is.null(scenario) || (scenario == "")) return()
-  cc <- paste(prefix, "country", sep = "_")
-  countries <- get_countries_ts(touchstone1, touchstone2, disease, group,
-                                scenario)
-  country <- update_dropdown_keep(session, cc, countries, input[[cc]])
-  update_country_ts(session, touchstone1, touchstone2, disease, group,
-                    scenario, country, input, prefix)
-}
-
-update_scenario_mg <- function(session, touchstone, disease, group1, group2,
-                               scenario, input, prefix = "bmg") {
-  if (is.null(scenario) || (scenario == "")) return()
-  if (is.null(group1) || is.null(group2) || (group1 == "") || (group2 == "")) return()
-  cc <- paste(prefix, "country", sep = "_")
-  countries <- get_countries_mg(touchstone, disease, group1, group2,
-                                scenario)
-  country <- update_dropdown_keep(session, cc, countries, input[[cc]])
-  update_country_mg(session, touchstone, disease, group1, group2,
-                    scenario, country, input, prefix)
-}
-
-update_country <- function(session, touchstone, disease, group, scenario, country,
-                           input, prefix = "b") {
-  if (is.null(country) || (country == "")) return()
-  o <- paste(prefix, "outcome", sep = "_")
-  outcomes <- get_outcomes(touchstone, disease, group, scenario, country)
-  update_dropdown_keep(session, o, outcomes, input[[o]])
-}
-
-update_country_ts <- function(session, touchstone1, touchstone2, disease, group,
-                              scenario, country, input, prefix = "bts") {
-  if (is.null(country) || (country == "")) return()
-  o <- paste(prefix, "outcome", sep = "_")
-  outcomes <- get_outcomes_ts(touchstone1, touchstone2, disease, group,
-                              scenario, country)
-  update_dropdown_keep(session, o, outcomes, input[[o]])
+  update_country(session, ts1, ts2, disease, g1, g2, s1, s2, country,
+                 input, prefix)
 }
 
-update_country_mg <- function(session, touchstone, disease, group1, group2,
-                              scenario, country, input, prefix = "bmg") {
+update_country <- function(session, ts1, ts2, disease, g1, g2, s1, s2, country,
+                           input, prefix) {
   if (is.null(country) || (country == "")) return()
-  if (is.null(group1) || is.null(group2) || (group1 == "") || (group2 == "")) return()
-  if (is.null(scenario) || (scenario == "")) return()
+  if (is.null(g1) || (g1 == "")) return()
   o <- paste(prefix, "outcome", sep = "_")
-  outcomes <- get_outcomes_mg(touchstone, disease, group1, group2,
-                              scenario, country)
+  outcomes <- get_outcomes(ts1, ts2, disease, g1, g2, s1, s2, country)
   update_dropdown_keep(session, o, outcomes, input[[o]])
 }

From 562259a1749429fb7e7bb739f1ff77f420d09160 Mon Sep 17 00:00:00 2001
From: Wes Hinsley <w.hinsley@imperial.ac.uk>
Date: Thu, 26 Mar 2026 20:39:46 +0000
Subject: [PATCH 11/29] Missing httr2 package

---
 DESCRIPTION | 1 +
 1 file changed, 1 insertion(+)

diff --git a/DESCRIPTION b/DESCRIPTION
index 1ab97ff..2a5b88e 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -20,6 +20,7 @@ Imports:
     cli,
     DBI,
     data.table,
+    httr2,
     jsonlite,
     lgr,
     dplyr,

From 805424940403199da6d9cff8d32dbc02a6ec324f Mon Sep 17 00:00:00 2001
From: Wes Hinsley <w.hinsley@imperial.ac.uk>
Date: Thu, 26 Mar 2026 20:50:07 +0000
Subject: [PATCH 12/29] Prevent autoplot

---
 inst/app/app.R | 31 ++++++++++++++++++++-----------
 1 file changed, 20 insertions(+), 11 deletions(-)

diff --git a/inst/app/app.R b/inst/app/app.R
index d1ed9e5..4f7a214 100644
--- a/inst/app/app.R
+++ b/inst/app/app.R
@@ -286,7 +286,16 @@ app_server <- function(input, output, session) {
 
       list(
         ages = ages,
-        opts = input[[sprintf("%s_options", prefix)]]
+        opts = input[[sprintf("%s_options", prefix)]],
+
+        touchstone = c(input[[it1]], if (!is.null(it2)) input[[it2]] else NULL),
+        group      = c(input[[ig1]], if (!is.null(ig2)) input[[ig2]] else NULL),
+        scenario   = c(input[[is1]], if (!is.null(is2)) input[[is2]] else NULL),
+
+        disease = input[[id]],
+        country = input[[ic]],
+        outcome = input[[sprintf("%s_outcome", prefix)]],
+        year    = input[[sprintf("%s_year", prefix)]]
       )
     })
 
@@ -298,19 +307,19 @@ app_server <- function(input, output, session) {
           pr <- plot_reactive()
           req(pr)
 
-          it <- if (n_touchstone == 1) it1 else if (gg == 1) it1 else it2
-          ig <- if (n_group == 1) ig1 else if (gg == 1) ig1 else ig2
+          it <- pr$touchstone[min(gg, length(pr$touchstone))]
+          ig <- pr$group[min(gg, length(pr$group))]
 
           stone_stochastic_graph(
             base = data_dir,
-            touchstone = input[[it]],
-            disease = input[[id]],
-            group = input[[ig]],
-            country = input[[ic]],
-            scenario = input[[is1]],
-            scenario2 = if (is.null(is2)) NULL else input[[is2]],
-            outcome = input[[sprintf("%s_outcome", prefix)]],
-            by_cohort = input[[sprintf("%s_year", prefix)]] == "Cohort",
+            touchstone = it,
+            disease = pr$disease,
+            group = ig,
+            country = pr$country,
+            scenario = pr$scenario[1],
+            scenario2 = if (length(pr$scenario) > 1) pr$scenario[2] else NULL,
+            outcome = pr$outcome,
+            by_cohort = pr$year == "Cohort",
             ages = pr$ages,
             log = "Log-Y" %in% pr$opts,
             include_median = "Median" %in% pr$opts,

From a9567615666f8f9a848d2a94392eeddd2e00d63d Mon Sep 17 00:00:00 2001
From: Wes Hinsley <w.hinsley@imperial.ac.uk>
Date: Thu, 26 Mar 2026 21:04:30 +0000
Subject: [PATCH 13/29] Parameterise data_dir

---
 R/stochastic_graphs.R      | 5 ++++-
 inst/app/app.R             | 2 --
 inst/app/utils.R           | 2 --
 man/stochastic_explorer.Rd | 2 +-
 4 files changed, 5 insertions(+), 6 deletions(-)

diff --git a/R/stochastic_graphs.R b/R/stochastic_graphs.R
index b6723c4..c70ab1a 100644
--- a/R/stochastic_graphs.R
+++ b/R/stochastic_graphs.R
@@ -160,6 +160,9 @@ prepare_central_data <- function(packit_id, packit_file,
 ##'
 ##' @export
 ##' @title Stochastic plot
-stochastic_explorer <- function() {
+##'
+stochastic_explorer <- function(
+  data_dir = "//wpia-hn2.hpc.dide.ic.ac.uk/vimc_stochastics") {
+  assign("data_dir", data_dir, envir = .GlobalEnv)
   shiny::runApp(system.file("app", package = "stoner"))
 }
diff --git a/inst/app/app.R b/inst/app/app.R
index 4f7a214..933819b 100644
--- a/inst/app/app.R
+++ b/inst/app/app.R
@@ -2,8 +2,6 @@ library(shiny)
 library(shinyjs)
 library(arrow)
 
-data_dir <- "//wpia-hn2.hpc.dide.ic.ac.uk/vimc_stochastics"
-
 source("utils.R")
 
 app_ui <- function() {
diff --git a/inst/app/utils.R b/inst/app/utils.R
index 1500413..daaa9cc 100644
--- a/inst/app/utils.R
+++ b/inst/app/utils.R
@@ -7,8 +7,6 @@
 # Functions ending in _mg will find things common
 # between to modelling groups, for comparisons.
 
-data_dir <- "//wpia-hn2.hpc.dide.ic.ac.uk/vimc_stochastics"
-
 get_touchstones <- function() {
   sort(unique(basename(list.dirs(data_dir, recursive = FALSE))))
 }
diff --git a/man/stochastic_explorer.Rd b/man/stochastic_explorer.Rd
index 5af6cd7..5649eb6 100644
--- a/man/stochastic_explorer.Rd
+++ b/man/stochastic_explorer.Rd
@@ -4,7 +4,7 @@
 \alias{stochastic_explorer}
 \title{Stochastic plot}
 \usage{
-stochastic_explorer()
+stochastic_explorer(data_dir = "//wpia-hn2.hpc.dide.ic.ac.uk/vimc_stochastics")
 }
 \description{
 Launch a Shiny app to allow interactive plotting of

From a2e08306fa75593cf9730660acea1fb9aaa45c6b Mon Sep 17 00:00:00 2001
From: Wes Hinsley <w.hinsley@imperial.ac.uk>
Date: Thu, 26 Mar 2026 21:12:27 +0000
Subject: [PATCH 14/29] Import warings

---
 NAMESPACE             | 2 ++
 R/stochastic_graphs.R | 3 ++-
 2 files changed, 4 insertions(+), 1 deletion(-)

diff --git a/NAMESPACE b/NAMESPACE
index e9fe275..853aae6 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -23,7 +23,9 @@ import(readr)
 importFrom(arrow,write_parquet)
 importFrom(data.table,as.data.table)
 importFrom(data.table,rbindlist)
+importFrom(grDevices,recordPlot)
 importFrom(graphics,lines)
+importFrom(graphics,par)
 importFrom(magrittr,"%>%")
 importFrom(rlang,":=")
 importFrom(stats,median)
diff --git a/R/stochastic_graphs.R b/R/stochastic_graphs.R
index c70ab1a..7b5c896 100644
--- a/R/stochastic_graphs.R
+++ b/R/stochastic_graphs.R
@@ -6,8 +6,9 @@
 ##' @import dplyr
 ##' @importFrom rlang :=
 ##' @import arrow
-##' @importFrom graphics lines
+##' @importFrom graphics lines par
 ##' @importFrom stats quantile median
+##' @importFrom grDevices recordPlot
 ##' @param base The folder in which the standardised stochastic files are found.
 ##' @param touchstone The touchstone name (for the graph title)
 ##' @param disease The disease, used for building the filename and graph title.

From c52d59b698df448f7f9ff169d017320088e5b6a6 Mon Sep 17 00:00:00 2001
From: Wes Hinsley <w.hinsley@imperial.ac.uk>
Date: Fri, 27 Mar 2026 11:13:45 +0000
Subject: [PATCH 15/29] Fix package warnings

---
 NAMESPACE                           | 1 +
 R/packit.R                          | 1 -
 R/stochastic_files.R                | 2 ++
 R/stochastic_graphs.R               | 9 ++++++---
 man/stochastic_explorer.Rd          | 5 +++++
 man/stone_stochastic_standardise.Rd | 3 +++
 6 files changed, 17 insertions(+), 4 deletions(-)

diff --git a/NAMESPACE b/NAMESPACE
index 853aae6..0fa8961 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -28,6 +28,7 @@ importFrom(graphics,lines)
 importFrom(graphics,par)
 importFrom(magrittr,"%>%")
 importFrom(rlang,":=")
+importFrom(shiny,runApp)
 importFrom(stats,median)
 importFrom(stats,quantile)
 importFrom(testthat,expect_equal)
diff --git a/R/packit.R b/R/packit.R
index 03177cd..c35afb4 100644
--- a/R/packit.R
+++ b/R/packit.R
@@ -12,7 +12,6 @@
 ##' @param server By default, the URL to the packit API on Montagu,
 ##' but this can be set to other packit API's if we want.
 ##' @returns The filename of the temporary file which has been downloaded.
-
 fetch_packit <- function(packet_id, filename,
   server = "https://montagu.vaccineimpact.org/packit/api/") {
 
diff --git a/R/stochastic_files.R b/R/stochastic_files.R
index 0ed805c..4322a2d 100644
--- a/R/stochastic_files.R
+++ b/R/stochastic_files.R
@@ -43,6 +43,8 @@
 ##' process.
 ##' @param allow_missing_yll yll was introduced in 2023? This flag allows
 ##' it to be missing for processing older stochastics.
+##' @param allow_missing_dalys Some early groups did not provide dalys; this
+##' flag allows dalys to be skipped.
 ##' @param allow_missing_indexes In some early runs, different groups
 ##' provided different numbers of files for each scenario, because some
 ##' countries did not implement particular coverage campaigns. This
diff --git a/R/stochastic_graphs.R b/R/stochastic_graphs.R
index 7b5c896..840609f 100644
--- a/R/stochastic_graphs.R
+++ b/R/stochastic_graphs.R
@@ -140,10 +140,10 @@ prepare_central_data <- function(packit_id, packit_file,
   central <- central[central$scenario == scenario, ]
   central <- central[central$burden_outcome == outcome, ]
   if (!is.null(ages)) {
-    central <- central[d$age %in% ages, ]
+    central <- central[central$age %in% ages, ]
   }
   if (by_cohort) {
-    central$year <- central$year - d$age
+    central$year <- central$year - central$age
   }
   names(central)[names(central) == "value"] <- outcome
   central <- central %>% group_by(.data$year) %>%
@@ -160,8 +160,11 @@ prepare_central_data <- function(packit_id, packit_file,
 ##' comparisons between modelling groups.
 ##'
 ##' @export
+##' @importFrom shiny runApp
 ##' @title Stochastic plot
-##'
+##' @param data_dir The location of the standardised stochastic folder
+##' hierarchy; this can be a local path, a fully-qualified network path on
+##' windows, or a mount point on linux or Mac.
 stochastic_explorer <- function(
   data_dir = "//wpia-hn2.hpc.dide.ic.ac.uk/vimc_stochastics") {
   assign("data_dir", data_dir, envir = .GlobalEnv)
diff --git a/man/stochastic_explorer.Rd b/man/stochastic_explorer.Rd
index 5649eb6..f3683ca 100644
--- a/man/stochastic_explorer.Rd
+++ b/man/stochastic_explorer.Rd
@@ -6,6 +6,11 @@
 \usage{
 stochastic_explorer(data_dir = "//wpia-hn2.hpc.dide.ic.ac.uk/vimc_stochastics")
 }
+\arguments{
+\item{data_dir}{The location of the standardised stochastic folder
+hierarchy; this can be a local path, a fully-qualified network path on
+windows, or a mount point on linux or Mac.}
+}
 \description{
 Launch a Shiny app to allow interactive plotting of
 standardised stochastic data, burden estimates,
diff --git a/man/stone_stochastic_standardise.Rd b/man/stone_stochastic_standardise.Rd
index dcf42ae..dadbaa6 100644
--- a/man/stone_stochastic_standardise.Rd
+++ b/man/stone_stochastic_standardise.Rd
@@ -68,6 +68,9 @@ process.}
 \item{allow_missing_yll}{yll was introduced in 2023? This flag allows
 it to be missing for processing older stochastics.}
 
+\item{allow_missing_dalys}{Some early groups did not provide dalys; this
+flag allows dalys to be skipped.}
+
 \item{allow_missing_indexes}{In some early runs, different groups
 provided different numbers of files for each scenario, because some
 countries did not implement particular coverage campaigns. This

From d6e78be60851a9a261d5832d286d5f88a33d5a9b Mon Sep 17 00:00:00 2001
From: Wes Hinsley <w.hinsley@imperial.ac.uk>
Date: Fri, 27 Mar 2026 11:21:11 +0000
Subject: [PATCH 16/29] Missing pkg

---
 DESCRIPTION | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/DESCRIPTION b/DESCRIPTION
index 2a5b88e..44c6151 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -28,6 +28,8 @@ Imports:
     prettyunits,
     readr,
     rlang,
+    shiny,
+    shinyjs,
     testthat,
     utils,
     withr

From 186ccf47b4a104f70198947ce9421908b82fb5b7 Mon Sep 17 00:00:00 2001
From: Wes Hinsley <w.hinsley@imperial.ac.uk>
Date: Fri, 27 Mar 2026 11:29:31 +0000
Subject: [PATCH 17/29] Packit/Packet typo

---
 R/packit.R          | 6 +++---
 man/fetch_packit.Rd | 8 ++++----
 2 files changed, 7 insertions(+), 7 deletions(-)

diff --git a/R/packit.R b/R/packit.R
index c35afb4..356cdbc 100644
--- a/R/packit.R
+++ b/R/packit.R
@@ -4,11 +4,11 @@
 ##' file from packet (ie, from the Montagu Reporting Portal).
 ##'
 ##' @export
-##' @title Fetch packit
+##' @title Fetch packet from a packit server
 ##' @import httr2
 ##' @importFrom tools file_ext
-##' @param packit_id The id of the packit containing the artefact.
-##' @param filename The filename of the file within the packit.
+##' @param packet_id The id of the packet containing the artefact.
+##' @param filename The filename of the file within the packet.
 ##' @param server By default, the URL to the packit API on Montagu,
 ##' but this can be set to other packit API's if we want.
 ##' @returns The filename of the temporary file which has been downloaded.
diff --git a/man/fetch_packit.Rd b/man/fetch_packit.Rd
index 3326eca..e77e94f 100644
--- a/man/fetch_packit.Rd
+++ b/man/fetch_packit.Rd
@@ -2,7 +2,7 @@
 % Please edit documentation in R/packit.R
 \name{fetch_packit}
 \alias{fetch_packit}
-\title{Fetch packit}
+\title{Fetch packet from a packit server}
 \usage{
 fetch_packit(
   packet_id,
@@ -11,12 +11,12 @@ fetch_packit(
 )
 }
 \arguments{
-\item{filename}{The filename of the file within the packit.}
+\item{packet_id}{The id of the packet containing the artefact.}
+
+\item{filename}{The filename of the file within the packet.}
 
 \item{server}{By default, the URL to the packit API on Montagu,
 but this can be set to other packit API's if we want.}
-
-\item{packit_id}{The id of the packit containing the artefact.}
 }
 \value{
 The filename of the temporary file which has been downloaded.

From 41b2f12165687330f5037c0722872f32fc63afde Mon Sep 17 00:00:00 2001
From: Wes Hinsley <w.hinsley@imperial.ac.uk>
Date: Fri, 27 Mar 2026 14:12:03 +0000
Subject: [PATCH 18/29] stoner:: on graph call

---
 inst/app/app.R | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/inst/app/app.R b/inst/app/app.R
index 933819b..612c305 100644
--- a/inst/app/app.R
+++ b/inst/app/app.R
@@ -307,8 +307,9 @@ app_server <- function(input, output, session) {
 
           it <- pr$touchstone[min(gg, length(pr$touchstone))]
           ig <- pr$group[min(gg, length(pr$group))]
+          browser()
 
-          stone_stochastic_graph(
+          stoner::stone_stochastic_graph(
             base = data_dir,
             touchstone = it,
             disease = pr$disease,

From 5815411f7623c298255af615fcf58150159b5362 Mon Sep 17 00:00:00 2001
From: Wes Hinsley <w.hinsley@imperial.ac.uk>
Date: Fri, 27 Mar 2026 14:36:16 +0000
Subject: [PATCH 19/29] Add spinny icon

---
 inst/app/app.R | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/inst/app/app.R b/inst/app/app.R
index 612c305..bef6258 100644
--- a/inst/app/app.R
+++ b/inst/app/app.R
@@ -65,6 +65,9 @@ app_ui <- function() {
 
   fluidPage(
     tags$style(HTML("
+      .shiny-busy {
+        cursor: wait !important;
+      }
       .selectize-input {
         white-space: nowrap;
         overflow: hidden;
@@ -307,7 +310,6 @@ app_server <- function(input, output, session) {
 
           it <- pr$touchstone[min(gg, length(pr$touchstone))]
           ig <- pr$group[min(gg, length(pr$group))]
-          browser()
 
           stoner::stone_stochastic_graph(
             base = data_dir,

From ef0d54f473ed0e4b8e1c8a7cb0fe05976bcf876f Mon Sep 17 00:00:00 2001
From: Wes Hinsley <w.hinsley@imperial.ac.uk>
Date: Fri, 27 Mar 2026 15:04:41 +0000
Subject: [PATCH 20/29] Typo

---
 R/stochastic_files.R                | 4 ++--
 man/stone_stochastic_standardise.Rd | 4 ++--
 2 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/R/stochastic_files.R b/R/stochastic_files.R
index 4322a2d..1ed001b 100644
--- a/R/stochastic_files.R
+++ b/R/stochastic_files.R
@@ -31,10 +31,10 @@
 ##' these to the simpler names. Processing Rubella stochastic files without
 ##' this set to TRUE will fail - so while we should always do this, keeping
 ##' the parameter makes it more clear in the code what we're doing and why.
-##' @param hepb2019_fix In 2019, HepB deaths and cases were subdivided into a
+##' @param hepb2019_fix In 2019 (and 2017), HepB deaths and cases were subdivided into
 ##' number of different causes. This flag combines those into the single
 ##' appropriate burden outcome.
-##' @param hib2019_fix In 2019, HepB deaths and cases were subdivided into a
+##' @param hib2019_fix In 2019 (and 2017), Hib deaths and cases were subdivided into
 ##' number of different causes. This flag combines those into the single
 ##' appropriate burden outcome.
 ##' @param missing_run_id_fix Some groups in the past have omitted run_id
diff --git a/man/stone_stochastic_standardise.Rd b/man/stone_stochastic_standardise.Rd
index dadbaa6..d010e20 100644
--- a/man/stone_stochastic_standardise.Rd
+++ b/man/stone_stochastic_standardise.Rd
@@ -52,11 +52,11 @@ these to the simpler names. Processing Rubella stochastic files without
 this set to TRUE will fail - so while we should always do this, keeping
 the parameter makes it more clear in the code what we're doing and why.}
 
-\item{hepb2019_fix}{In 2019, HepB deaths and cases were subdivided into a
+\item{hepb2019_fix}{In 2019 (and 2017), HepB deaths and cases were subdivided into
 number of different causes. This flag combines those into the single
 appropriate burden outcome.}
 
-\item{hib2019_fix}{In 2019, HepB deaths and cases were subdivided into a
+\item{hib2019_fix}{In 2019 (and 2017), Hib deaths and cases were subdivided into
 number of different causes. This flag combines those into the single
 appropriate burden outcome.}
 

From e34909d692cb6b50b290d8393ccb943803731847 Mon Sep 17 00:00:00 2001
From: Wes Hinsley <w.hinsley@imperial.ac.uk>
Date: Sat, 28 Mar 2026 19:20:14 +0000
Subject: [PATCH 21/29] Warn on file problem on startup

---
 R/stochastic_graphs.R | 12 ++++++++++++
 1 file changed, 12 insertions(+)

diff --git a/R/stochastic_graphs.R b/R/stochastic_graphs.R
index 840609f..c617ded 100644
--- a/R/stochastic_graphs.R
+++ b/R/stochastic_graphs.R
@@ -167,6 +167,18 @@ prepare_central_data <- function(packit_id, packit_file,
 ##' windows, or a mount point on linux or Mac.
 stochastic_explorer <- function(
   data_dir = "//wpia-hn2.hpc.dide.ic.ac.uk/vimc_stochastics") {
+  if (!dir.exists(file.path(data_dir))) {
+    cli::cli_abort(c(
+      "x" = "Cannot access the path/mount: {.path {data_dir}}",
+      "i" = "Please check you can see this path normally. If not:",
+      "*" = "You need ZScaler on and connected if you're not within DIDE",
+      "*" = "(But strongly advise remote desktop into DIDE - large files)",
+      "*" = "On linux, ensure the mount is correctly set up.",
+      "*" = "Check with DIDE IT that you have access to VIMC files",
+      "*" = "Check your general internet access."
+    ))
+  }
+
   assign("data_dir", data_dir, envir = .GlobalEnv)
   shiny::runApp(system.file("app", package = "stoner"))
 }

From d6ad5bc1f99ee79cd538be794d52cc3eb3a0fc2b Mon Sep 17 00:00:00 2001
From: Wes Hinsley <w.hinsley@imperial.ac.uk>
Date: Mon, 30 Mar 2026 14:21:00 +0100
Subject: [PATCH 22/29] Fix missing IC Hepb deaths outcome and test

---
 R/stochastic_files.R                   |   2 +-
 tests/testthat/test_stochastic_files.R | 130 +++++++++++++++++++++++++
 2 files changed, 131 insertions(+), 1 deletion(-)

diff --git a/R/stochastic_files.R b/R/stochastic_files.R
index 1ed001b..a185813 100644
--- a/R/stochastic_files.R
+++ b/R/stochastic_files.R
@@ -122,7 +122,7 @@ stone_stochastic_standardise <- function(
           ic_deaths <- c("hepb_deaths_acute", "hepb_deaths_comp_cirrh",
                          "hepb_deaths_dec_cirrh", "hepb_deaths_hcc")
 
-          for (i in unique(c(cda_deaths, li_deaths))) {
+          for (i in unique(c(cda_deaths, li_deaths, ic_deaths))) {
             if (i %in% names(d)) {
               message(sprintf("Including %s in deaths", i))
               d$deaths <- d$deaths + d[[i]]
diff --git a/tests/testthat/test_stochastic_files.R b/tests/testthat/test_stochastic_files.R
index d1a668a..b603140 100644
--- a/tests/testthat/test_stochastic_files.R
+++ b/tests/testthat/test_stochastic_files.R
@@ -147,3 +147,133 @@ test_that("Create central works", {
 
   expect_true(p1$cases[p1$country == "NPL"] == 2 * p1$cases[p1$country == "LAP"])
 })
+
+test_that("Rubella fix works", {
+  df <- fake_data()
+  df$rubella_deaths_congenital <- df$deaths
+  df$rubella_cases_congenital <- df$cases
+  df$deaths <- NULL
+  df$cases <- NULL
+  df$rubella_infections <- sample(nrow(df))
+
+  tmpin <- tempdir()
+  tmpout <- tempdir()
+  tmpfile <- tempfile(tmpdir = tmpin)
+  write.csv(df, paste0(tmpfile, "_opt"), row.names = FALSE)
+
+  stone_stochastic_standardise(
+    group = "north_pole_rub", in_path = tmpin, out_path = tmpout,
+    scenarios = "opt",
+    files = paste0(basename(tmpfile), "_opt"),
+  )
+
+  files <- list.files(path = tmpout)
+  expect_true("north_pole_rub_opt_LAP.pq" %in% files)
+  pq <- arrow::read_parquet(file.path(tmpout, "north_pole_rub_opt_LAP.pq"))
+  expect_true("cases" %in% names(pq))
+  expect_true("deaths" %in% names(pq))
+  expect_false("rubella_deaths_congenital" %in% names(pq))
+  expect_false("rubella_cases_congenital" %in% names(pq))
+  expect_false("rubella_infections" %in% names(pq))
+  expect_true(all(pq$cases == df$rubella_cases_congenital))
+  expect_true(all(pq$deaths == df$rubella_deaths_congenital))
+})
+
+test_that("Hib/PCV fix works", {
+  df <- fake_data()
+  df$cases_pneumo <- df$cases + 1
+  df$cases_men <- df$cases + 2
+  df$cases <- NULL
+  df$deaths_pneumo <- df$deaths + 3
+  df$deaths_men <- df$deaths + 4
+  df$deaths <- NULL
+
+  tmpin <- tempdir()
+  tmpout <- tempdir()
+  tmpfile <- tempfile(tmpdir = tmpin)
+  write.csv(df, paste0(tmpfile, "_opt"), row.names = FALSE)
+
+  stone_stochastic_standardise(
+    group = "north_pole_hib", in_path = tmpin, out_path = tmpout,
+    scenarios = "opt",
+    files = paste0(basename(tmpfile), "_opt"),
+  )
+
+  files <- list.files(path = tmpout)
+  expect_true("north_pole_hib_opt_LAP.pq" %in% files)
+  pq <- arrow::read_parquet(file.path(tmpout, "north_pole_hib_opt_LAP.pq"))
+  expect_true("cases" %in% names(pq))
+  expect_true("deaths" %in% names(pq))
+  expect_false("cases_men" %in% names(pq))
+  expect_false("cases_pneumo" %in% names(pq))
+  expect_false("deaths_men" %in% names(pq))
+  expect_false("deaths_pneumo" %in% names(pq))
+
+  expect_true(all(pq$cases == df$cases_pneumo + df$cases_men))
+  expect_true(all(pq$deaths == df$deaths_pneumo + df$deaths_men))
+
+})
+
+test_that("HepB fix works", {
+  df <- fake_data()
+  df$hepb_cases_acute_severe <- df$cases + 1
+  df$hepb_cases_dec_cirrh <- df$cases + 2
+  df$hepb_cases_hcc <- df$cases + 3
+  df$hepb_cases_acute_symp <- df$cases + 4
+  df$hepb_cases_fulminant <- df$cases + 5
+  df$hepb_cases_chronic <- df$cases + 6
+  df$hepb_chronic_symptomatic_in_acute_phase <- df$cases + 7
+  df$hepb_cases_comp_cirrh <- df$cases + 8
+  df$hepb_cases_hcc_no_cirrh <- df$cases + 9
+  df$cases <- NULL
+  df$hepb_deaths_acute <- df$deaths + 1
+  df$hepb_deaths_dec_cirrh <- df$deaths + 2
+  df$hepb_deaths_hcc <- df$deaths + 3
+  df$hepb_deaths_total_cirrh <- df$deaths + 4
+  df$hepb_deaths_comp_cirrh <- df$deaths + 5
+  df$deaths <- NULL
+
+  tmpin <- tempdir()
+  tmpout <- tempdir()
+  tmpfile <- tempfile(tmpdir = tmpin)
+  write.csv(df, paste0(tmpfile, "_opt"), row.names = FALSE)
+
+  stone_stochastic_standardise(
+    group = "north_pole_hepb", in_path = tmpin, out_path = tmpout,
+    scenarios = "opt",
+    files = paste0(basename(tmpfile), "_opt"),
+  )
+
+  files <- list.files(path = tmpout)
+  expect_true("north_pole_hepb_opt_LAP.pq" %in% files)
+  pq <- arrow::read_parquet(file.path(tmpout, "north_pole_hepb_opt_LAP.pq"))
+
+  expect_true("cases" %in% names(pq))
+  expect_true("deaths" %in% names(pq))
+  expect_false("hepb_cases_acute_severe" %in% names(pq))
+  expect_false("hepb_cases_dec_cirrh" %in% names(pq))
+  expect_false("hepb_cases_hcc" %in% names(pq))
+  expect_false("hepb_cases_acute_symp" %in% names(pq))
+  expect_false("hepb_cases_fulminant" %in% names(pq))
+  expect_false("hepb_cases_chronic" %in% names(pq))
+  expect_false("hepb_chronic_symptomatic_in_acute_phase" %in% names(pq))
+  expect_false("hepb_cases_comp_cirrh" %in% names(pq))
+  expect_false("hepb_cases_comp_hcc_no_cirrh" %in% names(pq))
+
+  expect_false("hepb_deaths_acute" %in% names(pq))
+  expect_false("hepb_deaths_dec_cirrh" %in% names(pq))
+  expect_false("hepb_deaths_hcc" %in% names(pq))
+  expect_false("hepb_deaths_total_cirrh" %in% names(pq))
+  expect_false("hepb_deaths_comp_cirrh" %in% names(pq))
+
+  expect_true(all(pq$cases == df$hepb_cases_acute_severe + df$hepb_cases_dec_cirrh +
+                          df$hepb_Cases_hcc + df$hepb_cases_acute_symp +
+                          df$hepb_cases_fulminant + df$hepb_cases_chronic +
+                          df$hepb_chronic_symptomatic_in_acute_phase +
+                          df$hepb_cases_comp_cirrh + df$hepb_cases_comp_hcc_no_cirrh))
+
+  expect_true(all(pq$deaths == df$hepb_deaths_acute + df$hepb_deaths_dec_cirrh +
+                           df$hep_deaths_hcc + df$hepb_deaths_total_cirrh +
+                           df$hepb_deaths_comp_cirrh))
+
+})

From 83fe5d8064606777048187b44ca599984a4d40a7 Mon Sep 17 00:00:00 2001
From: Wes Hinsley <w.hinsley@imperial.ac.uk>
Date: Mon, 30 Mar 2026 15:00:28 +0100
Subject: [PATCH 23/29] More tests

---
 R/stochastic_graphs.R                   |  2 +-
 tests/testthat/test_stochastic_files.R  | 59 +++++++++++++++++++++++--
 tests/testthat/test_stochastic_graphs.R | 29 ++++++++++++
 3 files changed, 86 insertions(+), 4 deletions(-)

diff --git a/R/stochastic_graphs.R b/R/stochastic_graphs.R
index c617ded..76fbf07 100644
--- a/R/stochastic_graphs.R
+++ b/R/stochastic_graphs.R
@@ -180,5 +180,5 @@ stochastic_explorer <- function(
   }
 
   assign("data_dir", data_dir, envir = .GlobalEnv)
-  shiny::runApp(system.file("app", package = "stoner"))
+  runApp(system.file("app", package = "stoner"))
 }
diff --git a/tests/testthat/test_stochastic_files.R b/tests/testthat/test_stochastic_files.R
index b603140..ac19da8 100644
--- a/tests/testthat/test_stochastic_files.R
+++ b/tests/testthat/test_stochastic_files.R
@@ -164,7 +164,7 @@ test_that("Rubella fix works", {
   stone_stochastic_standardise(
     group = "north_pole_rub", in_path = tmpin, out_path = tmpout,
     scenarios = "opt",
-    files = paste0(basename(tmpfile), "_opt"),
+    files = paste0(basename(tmpfile), "_opt")
   )
 
   files <- list.files(path = tmpout)
@@ -196,7 +196,7 @@ test_that("Hib/PCV fix works", {
   stone_stochastic_standardise(
     group = "north_pole_hib", in_path = tmpin, out_path = tmpout,
     scenarios = "opt",
-    files = paste0(basename(tmpfile), "_opt"),
+    files = paste0(basename(tmpfile), "_opt")
   )
 
   files <- list.files(path = tmpout)
@@ -241,7 +241,7 @@ test_that("HepB fix works", {
   stone_stochastic_standardise(
     group = "north_pole_hepb", in_path = tmpin, out_path = tmpout,
     scenarios = "opt",
-    files = paste0(basename(tmpfile), "_opt"),
+    files = paste0(basename(tmpfile), "_opt")
   )
 
   files <- list.files(path = tmpout)
@@ -275,5 +275,58 @@ test_that("HepB fix works", {
   expect_true(all(pq$deaths == df$hepb_deaths_acute + df$hepb_deaths_dec_cirrh +
                            df$hep_deaths_hcc + df$hepb_deaths_total_cirrh +
                            df$hepb_deaths_comp_cirrh))
+})
+
+test_that("Missing yll/dalys works", {
+  df <- fake_data()
+  df$dalys <- NULL
+  df$yll <- NULL
+  tmpin <- tempdir()
+  tmpout <- tempdir()
+  tmpfile <- tempfile(tmpdir = tmpin)
+  write.csv(df, paste0(tmpfile, "_opt"), row.names = FALSE)
+
+  stone_stochastic_standardise(
+    group = "north_pole_flu", in_path = tmpin, out_path = tmpout,
+    scenarios = "opt",
+    files = paste0(basename(tmpfile), "_opt")
+  )
+  files <- list.files(path = tmpout)
+  expect_true("north_pole_flu_opt_LAP.pq" %in% files)
+  pq <- arrow::read_parquet(file.path(tmpout, "north_pole_flu_opt_LAP.pq"))
+  expect_false("dalys" %in% names(pq))
+  expect_false("yll" %in% names(pq))
+})
 
+test_that("Different file count per scenario is handled", {
+  fake <- fake_data()
+  tmpin <- tempdir()
+  tmpout <- tempdir()
+  tmpfile <- tempfile(tmpdir = tmpin)
+
+  fake$country <- "POL"
+  write.csv(fake, paste0(tmpfile, "_optimistic_1"), row.names = FALSE)
+  write.csv(fake, paste0(tmpfile, "_fatalistic_1"), row.names = FALSE)
+  fake$country <- "NOR"
+  write.csv(fake, paste0(tmpfile, "_optimistic_2"), row.names = FALSE)
+  write.csv(fake, paste0(tmpfile, "_fatalistic_2"), row.names = FALSE)
+  fake$country <- "LAP"
+  write.csv(fake, paste0(tmpfile, "_fatalistic_3"), row.names = FALSE)
+
+  stone_stochastic_standardise(
+    group = "north_pole_lurgy",
+    in_path = tmpin,
+    out_path = tmpout,
+    scenarios = c("optimistic", "fatalistic"),
+    files = paste0(basename(tmpfile), "_:scenario_:index"),
+    index = 1:3,
+    allow_missing_indexes = TRUE
+  )
+  files <- list.files(path = tmpout)
+  expect_true("north_pole_lurgy_optimistic_POL.pq" %in% files)
+  expect_true("north_pole_lurgy_optimistic_NOR.pq" %in% files)
+  expect_false("north_pole_lurgy_optimistic_LAP.pq" %in% files)
+  expect_true("north_pole_lurgy_fatalistic_POL.pq" %in% files)
+  expect_true("north_pole_lurgy_fatalistic_NOR.pq" %in% files)
+  expect_true("north_pole_lurgy_fatalistic_LAP.pq" %in% files)
 })
diff --git a/tests/testthat/test_stochastic_graphs.R b/tests/testthat/test_stochastic_graphs.R
index 243d3d3..8da222a 100644
--- a/tests/testthat/test_stochastic_graphs.R
+++ b/tests/testthat/test_stochastic_graphs.R
@@ -65,5 +65,34 @@ test_that("stochastic_graph data transforms", {
   expect_no_error(stone_stochastic_graph(
     base, touchstone, disease, group, country,
     scenario, "deaths", log = TRUE))
+})
+
+test_that("stochastic_explorer data_dir handling", {
+  expect_error(stochastic_explorer("//localhost/potato",
+                      "Cannot access the path/mount"))
+})
 
+test_that("Can launch shiny app", {
+  fake_path <- tempdir()
+  runApp_called <- FALSE
+  runApp_arg <- NULL
+
+  local_mocked_bindings(
+    runApp = function(app_dir) {
+      runApp_called <<- TRUE
+      runApp_arg <<- app_dir
+      invisible(NULL)
+    },
+    .env = environment(stochastic_explorer)
+  )
+
+  withr::with_envvar(c(), {
+    if (exists("data_dir", envir = .GlobalEnv)) {
+      rm(data_dir, envir = .GlobalEnv)
+    }
+    stochastic_explorer(data_dir = fake_path)
+    expect_true(exists("data_dir", envir = .GlobalEnv))
+    expect_equal(get("data_dir", envir = .GlobalEnv), fake_path)
+    expect_true(runApp_called)
+  })
 })

From 4d10598fd0ff539bd9986d84b401a3c86b79c275 Mon Sep 17 00:00:00 2001
From: Wes Hinsley <w.hinsley@imperial.ac.uk>
Date: Mon, 30 Mar 2026 15:27:43 +0100
Subject: [PATCH 24/29] Tested except for packit

---
 R/stochastic_graphs.R                   | 16 ++++++++++------
 tests/testthat/test_stochastic_graphs.R | 14 ++++++++++++--
 2 files changed, 22 insertions(+), 8 deletions(-)

diff --git a/R/stochastic_graphs.R b/R/stochastic_graphs.R
index 76fbf07..d471efb 100644
--- a/R/stochastic_graphs.R
+++ b/R/stochastic_graphs.R
@@ -1,3 +1,11 @@
+age_string <- function(ages) {
+  if (is.null(ages)) {
+    "all ages"
+  } else if (identical(unique(sort(ages)), as.numeric(min(ages):max(ages)))) {
+    sprintf("age %d..%d", min(ages), max(ages))
+  } else "selected ages"
+}
+
 ##' Draw a stochastic plot showing all the different runs, with the mean,
 ##' median, 5% and 95% quantiles shown.
 ##'
@@ -48,11 +56,7 @@ stone_stochastic_graph <- function(base, touchstone, disease, group, country,
 
   d <- prepare_graph_data(base, touchstone, disease, group, country,
                          scenario, outcome, ages, by_cohort)
-  if (is.null(ages)) {
-    age <- "all ages"
-  } else if (all(sort(ages) == min(ages):max(ages))) {
-    age <- sprintf("age %d..%d", min(ages), max(ages))
-  } else age <- "selected ages"
+  age <- age_string(ages)
 
   title <- sprintf("%s, %s, %s, %s\n%s, %s\n", touchstone, disease, group, age,
                    scenario, country)
@@ -167,7 +171,7 @@ prepare_central_data <- function(packit_id, packit_file,
 ##' windows, or a mount point on linux or Mac.
 stochastic_explorer <- function(
   data_dir = "//wpia-hn2.hpc.dide.ic.ac.uk/vimc_stochastics") {
-  if (!dir.exists(file.path(data_dir))) {
+  if (!dir.exists(data_dir)) {
     cli::cli_abort(c(
       "x" = "Cannot access the path/mount: {.path {data_dir}}",
       "i" = "Please check you can see this path normally. If not:",
diff --git a/tests/testthat/test_stochastic_graphs.R b/tests/testthat/test_stochastic_graphs.R
index 8da222a..192d909 100644
--- a/tests/testthat/test_stochastic_graphs.R
+++ b/tests/testthat/test_stochastic_graphs.R
@@ -65,11 +65,15 @@ test_that("stochastic_graph data transforms", {
   expect_no_error(stone_stochastic_graph(
     base, touchstone, disease, group, country,
     scenario, "deaths", log = TRUE))
+
+  expect_no_error(stone_stochastic_graph(
+    base, touchstone, disease, group, country,
+    scenario, "deaths", scenario2 = scenario))
 })
 
 test_that("stochastic_explorer data_dir handling", {
-  expect_error(stochastic_explorer("//localhost/potato",
-                      "Cannot access the path/mount"))
+  expect_error(stochastic_explorer(file.path(tempdir(), "potato", "salad")),
+                      "Cannot access the path/mount")
 })
 
 test_that("Can launch shiny app", {
@@ -96,3 +100,9 @@ test_that("Can launch shiny app", {
     expect_true(runApp_called)
   })
 })
+
+test_that("Age formats are reasonable", {
+  expect_equal(age_string(NULL), "all ages")
+  expect_equal(age_string(c(5,4,3,2,1,5,4,3,2,1)), "age 1..5")
+  expect_equal(age_string(c(2,4,6,8)), "selected ages")
+})

From 7ef223a8a3c9fd0818ac35c8912c984bb9c16e1e Mon Sep 17 00:00:00 2001
From: Wes Hinsley <w.hinsley@imperial.ac.uk>
Date: Mon, 30 Mar 2026 16:44:19 +0100
Subject: [PATCH 25/29] Test central creation from packit

---
 R/packit.R                              |  2 +-
 R/stochastic_graphs.R                   |  6 +--
 tests/testthat/test_stochastic_graphs.R | 61 +++++++++++++++++++++++++
 3 files changed, 65 insertions(+), 4 deletions(-)

diff --git a/R/packit.R b/R/packit.R
index 356cdbc..087b098 100644
--- a/R/packit.R
+++ b/R/packit.R
@@ -61,7 +61,7 @@ fetch_packit <- function(packet_id, filename,
 
   # And finally, we download the file in chunks.
 
-  out <- tempfile(fileext = tools::file_ext(filename))
+  out <- tempfile(fileext = paste0(".", tools::file_ext(filename)))
   outfile <- file(out, open = "wb")
   con <- httr2::req_perform_connection(req, blocking = TRUE)
   while (!httr2::resp_stream_is_complete(con)) {
diff --git a/R/stochastic_graphs.R b/R/stochastic_graphs.R
index d471efb..172369f 100644
--- a/R/stochastic_graphs.R
+++ b/R/stochastic_graphs.R
@@ -79,9 +79,9 @@ stone_stochastic_graph <- function(base, touchstone, disease, group, country,
   log <- if (log) "y" else ""
 
   if (!is.null(packit_id)) {
-    central <- prepare_central_data(packit_id, packit_file)
-
-
+    central <- prepare_central_data(packit_id, packit_file,
+                                    country, scenario, outcome, ages, by_cohort)
+    # To be continued...
   }
   par(mar = c(5, 4, 5, 2))
   plot(ylab = outcome_ylab, xlab = if (by_cohort) "Birth Cohort" else "year",
diff --git a/tests/testthat/test_stochastic_graphs.R b/tests/testthat/test_stochastic_graphs.R
index 192d909..8040dd4 100644
--- a/tests/testthat/test_stochastic_graphs.R
+++ b/tests/testthat/test_stochastic_graphs.R
@@ -69,6 +69,29 @@ test_that("stochastic_graph data transforms", {
   expect_no_error(stone_stochastic_graph(
     base, touchstone, disease, group, country,
     scenario, "deaths", scenario2 = scenario))
+
+  # Packit gets called if needed
+
+  called <- FALSE
+  called_args <- NULL
+
+  local_mocked_bindings(
+    prepare_central_data = function(id, file) {
+      called <<- TRUE
+      called_args <<- list(id = id, file = file)
+      "fake_result"
+    },
+    .env = environment(stone_stochastic_graph)
+  )
+
+  expect_no_error(stone_stochastic_graph(
+    base, touchstone, disease, group, country,
+    scenario, "deaths", packit_id = "123",
+    packit_file = "file.csv"))
+
+  expect_true(called)
+  expect_equal(called_args$id, "123")
+  expect_equal(called_args$file, "file.csv")
 })
 
 test_that("stochastic_explorer data_dir handling", {
@@ -106,3 +129,41 @@ test_that("Age formats are reasonable", {
   expect_equal(age_string(c(5,4,3,2,1,5,4,3,2,1)), "age 1..5")
   expect_equal(age_string(c(2,4,6,8)), "selected ages")
 })
+
+test_that("Parsing central from packit works", {
+  # Packit gets called if needed
+
+  fake <- data.frame(
+    scenario_type = "RSV-rout", scenario = "RSV-rout",
+    year = c(rep(2000, 4), rep(2001, 4), rep(2000, 4), rep(2001, 4)),
+    age = c(rep(0, 8), rep(1, 8)),
+    country = "RFP",
+    burden_outcome = rep(c("cases", "dalys", "deaths", "yll"), 2),
+    value = 1:16)
+
+  rds <- tempfile(fileext = ".rds")
+  saveRDS(fake, rds)
+
+  local_mocked_bindings(
+    fetch_packit = function(id, file) {
+      rds
+    },
+    .env = environment(prepare_central_data)
+  )
+
+  res <- prepare_central_data("123", "file.csv",
+    "RFP", "RSV-rout", "deaths", 0:5, TRUE)
+
+  # Data in for death is: (year, age, deaths)
+  # 2000, 0, 3
+  # 2000, 1, 11
+  # 2001, 0, 7
+  # 2001, 1, 15
+  # For cohort - this should become...
+  # 1999, 11   (2000 year 1, were born in 1999)
+  # 2000, 18   (2000 year 0, and 2001 year 1 born in 2000)
+  # 2001, 7    (2001 year 0)
+
+  expect_true(all.equal(res$year, c(1999, 2000, 2001)))
+  expect_true(all.equal(res$deaths, c(11, 18, 7)))
+})

From c5fb28b609d538f1a63c367e37311e104a0651a8 Mon Sep 17 00:00:00 2001
From: Wes Hinsley <w.hinsley@imperial.ac.uk>
Date: Mon, 30 Mar 2026 17:19:39 +0100
Subject: [PATCH 26/29] Fix mocks

---
 tests/testthat/test_stochastic_graphs.R | 28 +++++++------------------
 1 file changed, 8 insertions(+), 20 deletions(-)

diff --git a/tests/testthat/test_stochastic_graphs.R b/tests/testthat/test_stochastic_graphs.R
index 8040dd4..58a2735 100644
--- a/tests/testthat/test_stochastic_graphs.R
+++ b/tests/testthat/test_stochastic_graphs.R
@@ -72,26 +72,18 @@ test_that("stochastic_graph data transforms", {
 
   # Packit gets called if needed
 
-  called <- FALSE
-  called_args <- NULL
-
-  local_mocked_bindings(
-    prepare_central_data = function(id, file) {
-      called <<- TRUE
-      called_args <<- list(id = id, file = file)
-      "fake_result"
-    },
-    .env = environment(stone_stochastic_graph)
-  )
+  fake_result <- mockery::mock("fake_result")
+  mockery::stub(stone_stochastic_graph, "prepare_central_data", fake_result)
 
   expect_no_error(stone_stochastic_graph(
     base, touchstone, disease, group, country,
     scenario, "deaths", packit_id = "123",
     packit_file = "file.csv"))
 
-  expect_true(called)
-  expect_equal(called_args$id, "123")
-  expect_equal(called_args$file, "file.csv")
+  mockery::expect_called(fake_result, 1)
+  mockery::expect_args(fake_result, 1, "123", "file.csv", country,
+                       scenario, "deaths", NULL, FALSE)
+
 })
 
 test_that("stochastic_explorer data_dir handling", {
@@ -144,12 +136,8 @@ test_that("Parsing central from packit works", {
   rds <- tempfile(fileext = ".rds")
   saveRDS(fake, rds)
 
-  local_mocked_bindings(
-    fetch_packit = function(id, file) {
-      rds
-    },
-    .env = environment(prepare_central_data)
-  )
+  fetch_fake <- function(id, file) rds
+  mockery::stub(prepare_central_data, "fetch_packit", fetch_fake)
 
   res <- prepare_central_data("123", "file.csv",
     "RFP", "RSV-rout", "deaths", 0:5, TRUE)

From c5c0bf11f9903d398ed311ae2b5eef0e1fcf42f5 Mon Sep 17 00:00:00 2001
From: Wes Hinsley <w.hinsley@imperial.ac.uk>
Date: Mon, 30 Mar 2026 17:40:52 +0100
Subject: [PATCH 27/29] Ignore fetch_packit for covr

---
 .covrignore | 4 ++++
 1 file changed, 4 insertions(+)
 create mode 100644 .covrignore

diff --git a/.covrignore b/.covrignore
new file mode 100644
index 0000000..98b7977
--- /dev/null
+++ b/.covrignore
@@ -0,0 +1,4 @@
+# This does the API call to packit - not 
+# easy to test without mocking all the network.
+
+R/fetch_packit.R

From b40e3aadad99239fdff89a5e7e41f16a0e7b2e84 Mon Sep 17 00:00:00 2001
From: Wes Hinsley <w.hinsley@imperial.ac.uk>
Date: Mon, 30 Mar 2026 17:48:10 +0100
Subject: [PATCH 28/29] Use correct filename...

---
 .covrignore | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.covrignore b/.covrignore
index 98b7977..ddb6d3d 100644
--- a/.covrignore
+++ b/.covrignore
@@ -1,4 +1,4 @@
 # This does the API call to packit - not 
 # easy to test without mocking all the network.
 
-R/fetch_packit.R
+R/packit.R

From 18efc9bb72cdc6f2255f5bf4e42387f0dd042a7b Mon Sep 17 00:00:00 2001
From: Wes Hinsley <w.hinsley@imperial.ac.uk>
Date: Tue, 31 Mar 2026 12:11:55 +0100
Subject: [PATCH 29/29] Add mockery to suggests

---
 DESCRIPTION | 1 +
 1 file changed, 1 insertion(+)

diff --git a/DESCRIPTION b/DESCRIPTION
index 44c6151..8eea4ac 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -38,6 +38,7 @@ RoxygenNote: 7.3.3
 Roxygen: list(markdown = TRUE)
 Suggests:
     knitr,
+    mockery,
     rcmdcheck,
     rmarkdown,
     RPostgres