From f3af9a031fdaa10d0582ab58dc1eedc9e1fd920b Mon Sep 17 00:00:00 2001
From: Andrew Shattock <andrewjames.shattock@unibas.ch>
Date: Tue, 16 Apr 2024 17:01:30 +0200
Subject: [PATCH] Remove figures not in manuscript nor supplement

---
 coverage.R   |   11 -
 external.R   |    9 -
 history.R    |   10 -
 impact.R     |    3 -
 launch.R     |    2 +-
 options.R    |    1 -
 plotting.R   | 1457 +-------------------------------------------------
 regression.R |    3 -
 results.R    |   40 +-
 9 files changed, 30 insertions(+), 1506 deletions(-)

diff --git a/coverage.R b/coverage.R
index db774d2..0e3eef8 100644
--- a/coverage.R
+++ b/coverage.R
@@ -70,14 +70,8 @@ prepare_coverage = function() {
   # Plot total number of FVP over time
   plot_total_fvps()
   
-  # Plot coverage density by disease
-  plot_coverage()
-  
   # Coverage data density by age
   plot_coverage_age_density()
-  
-  # Missing coverage data by country
-  plot_missing_data()
 }
 
 # ---------------------------------------------------------
@@ -550,14 +544,9 @@ smooth_static_fvps = function(coverage_dt) {
     fill(source, .direction = "updown") %>%
     arrange(d_v_a_id, country, year, age)
   
-  # ---- Diagnostic plots ----
-  
   # Save table for diagnostic plots
   save_table(smooth_dt, "smoothed_fvps")
   
-  # Diagnostic plots for FVPs smoothing
-  plot_smooth_fvps()
-  
   return(smoothed_coverage_dt)
 }
 
diff --git a/external.R b/external.R
index b817862..b2a2f81 100644
--- a/external.R
+++ b/external.R
@@ -35,20 +35,11 @@ run_external = function() {
   
   # ---- Data visualisation plots ----
   
-  # Plot outcomes from each external model
-  plot_external_models()
-  
   # Plot total number of FVP over time
   plot_total_fvps()
   
-  # Plot coverage density by disease
-  plot_coverage()
-  
   # Coverage data density by age
   plot_coverage_age_density()
-  
-  # Missing coverage data by country
-  plot_missing_data()
 }
 
 # ---------------------------------------------------------
diff --git a/history.R b/history.R
index faa89cd..9ac8b5c 100644
--- a/history.R
+++ b/history.R
@@ -329,16 +329,6 @@ run_history = function(metric) {
     # Save results to file
     save_rds(yll_time_dt, "history", "burden_averted_yll") 
   }
-  
-  # ---- Plot diagnostic figures ----
-  
-  # NOTE: Most results figures created in final module
-  
-  # Plot inital impact ratios used to back project
-  plot_impact_fvps(metric, scope = "initial")
-  
-  # Plot comparison of EPI50 outcomes vs VIMC outcomes
-  plot_vimc_comparison()
 }
 
 # ---------------------------------------------------------
diff --git a/impact.R b/impact.R
index a7e390a..85905ee 100644
--- a/impact.R
+++ b/impact.R
@@ -23,9 +23,6 @@ run_impact = function(metric) {
   # Prepare impact-FVP data to fit to
   data_dt = get_impact_data(metric)
   
-  # Exploratory plots of data used to fit impact functions
-  plot_impact_data(metric)
-  
   # ---- Model fitting ----
   
   message(" > Evaluating impact functions")
diff --git a/launch.R b/launch.R
index 74020a0..aece9fc 100644
--- a/launch.R
+++ b/launch.R
@@ -17,7 +17,7 @@ source("dependencies.R")
 message("Running EPI50 pipeline")
 
 # Define modules to be run
-run_modules = 8
+run_modules = 1 : 8
 
 # Set global options (see options.R)
 o = set_options(run_module = run_modules)
diff --git a/options.R b/options.R
index 635b04c..28f1463 100644
--- a/options.R
+++ b/options.R
@@ -57,7 +57,6 @@ set_options = function(run_module = NA) {
 
   # Turn figures on or off
   o$plot_inputs     = TRUE
-  o$plot_external   = TRUE
   o$plot_static     = TRUE
   o$plot_imputation = TRUE
   o$plot_impact     = TRUE
diff --git a/plotting.R b/plotting.R
index b34ec67..a7fe9e6 100644
--- a/plotting.R
+++ b/plotting.R
@@ -407,186 +407,6 @@ plot_total_fvps = function() {
   save_fig(g, "S07")
 }
 
-# ---------------------------------------------------------
-# Plot smoothed FVP for static model pathogens
-# ---------------------------------------------------------
-plot_smooth_fvps = function() {
-  
-  message("  - Plotting smoothed FVPs (static model pathogens)")
-  
-  # ---- Plot 1: data vs smoothing ----
-  
-  # Smoothing data - group by country and age
-  smooth_dt = table("smoothed_fvps") %>%
-    mutate(country_age = paste1(country, age), 
-           fvps_smooth = fvps_smooth / 1e6, 
-           fvps        = fvps / 1e6) %>%
-    format_d_v_a_name() %>%
-    select(d_v_a_name, country, country_age, 
-           year, fvps, fvps_smooth)
-  
-  # Plot the data with associated smoothing
-  g1 = ggplot(smooth_dt) +
-    aes(x = year, 
-        colour = country, 
-        group  = country_age) +
-    # Plot data...
-    geom_point(
-      mapping = aes(y = fvps),
-      show.legend = FALSE) +
-    # Plot smoothing...
-    geom_line(
-      mapping = aes(y = fvps_smooth),
-      show.legend = FALSE) +
-    # Simple faceting...
-    facet_wrap(~d_v_a_name) + 
-    # Prettify x axis...
-    scale_x_continuous(
-      expand = expansion(mult = c(0, 0)), 
-      breaks = pretty_breaks()) +  
-    # Prettify y axis...
-    scale_y_continuous(
-      name   = "Fully vaccinated people (in millions)", 
-      labels = comma,
-      limits = c(0, NA),
-      expand = expansion(mult = c(0, 0.05)))
-  
-  # Prettify theme
-  g1 = g1 + theme_classic() + 
-    theme(axis.title.x  = element_blank(),
-          axis.title.y  = element_text(
-            size = 20, margin = margin(l = 10, r = 20)),
-          axis.text     = element_text(size = 10),
-          # axis.text.x   = element_text(hjust = 1, angle = 50), 
-          axis.line     = element_blank(),
-          strip.text    = element_text(size = 14),
-          strip.background = element_blank(), 
-          panel.border  = element_rect(
-            linewidth = 0.5, fill = NA),
-          panel.spacing = unit(1, "lines"))
-  
-  # ---- Plot 2: smoothing error ----
-  
-  # Total error for each vaccine
-  diagnostic_dt = smooth_dt %>%
-    # Summarise for each vaccine...
-    group_by(d_v_a_name) %>%
-    summarise(fvps = sum(fvps),
-              fvps_smooth = sum(fvps_smooth)) %>%
-    ungroup() %>%
-    # Take the absolute difference...
-    mutate(diff = fvps - fvps_smooth,
-           abs  = abs(diff)) %>%
-    # Calculate total error...
-    mutate(err = round(100 * abs / fvps, 2), 
-           err = paste0("Error: ", err, "%")) %>%
-    # Append error to vaccine description...
-    mutate(d_v_a_name = paste0(d_v_a_name, "\n", err)) %>%
-    # Set plotting order by abs diff...
-    arrange(abs) %>%
-    mutate(d_v_a_name = fct_inorder(d_v_a_name)) %>%
-    as.data.table()
-  
-  # Plot total smoothing errors
-  g2 = ggplot(diagnostic_dt) +
-    aes(x = d_v_a_name, 
-        y = diff, 
-        fill = abs) +
-    geom_col(show.legend = FALSE) +
-    coord_flip() + 
-    # Prettify x axis (noting coord_flip)...
-    scale_y_continuous(
-      name   = "Total difference in FVPs after smoothing (in millions)", 
-      labels = comma)
-  
-  # Prettify theme
-  g2 = g2 + theme_classic() + 
-    theme(axis.text     = element_text(size = 12),
-          axis.title.x  = element_text(
-            size = 20, margin = margin(b = 10, t = 20)),
-          axis.title.y  = element_blank(),
-          axis.line     = element_blank(),
-          panel.border  = element_rect(
-            linewidth = 0.5, fill = NA))
-  
-  # ---- Save diagnostic plots ----
-  
-  # Save figures to file
-  save_fig(g1, "X03")
-  save_fig(g2, "X04")
-}
-
-# ---------------------------------------------------------
-# Plot coverage density by disease
-# ---------------------------------------------------------
-plot_coverage = function() {
-  
-  # All coverage values by disease, vaccine, and age class
-  all_coverage_dt = table("coverage") %>%
-    # Classify by age group...
-    mutate(age_group = ifelse(age == 0, "infant", "other")) %>%
-    # Classify by disease...
-    left_join(y  = table("d_v_a"), 
-              by = "d_v_a_id") %>%
-    # Tidy up...
-    select(disease, vaccine, age_group, coverage) %>%
-    arrange(disease, vaccine, age_group)
-  
-  # Produce plot for both disease and vaccine
-  for (d_v in c("disease", "vaccine")) {
-    
-    # Select appropriate column: disease or vaccine
-    plot_dt = all_coverage_dt %>%
-      select(d_v = !!d_v, age_group, coverage)
-    
-    # Plot coverage value density
-    g = ggplot(plot_dt) + 
-      aes(x = coverage, 
-          y = after_stat(scaled), 
-          colour = age_group,
-          fill   = age_group) + 
-      geom_density(alpha = 0.2) + 
-      facet_wrap(~d_v) + 
-      # Prettify x axis...
-      scale_x_continuous(
-        name   = "Coverage",
-        labels = percent,
-        limits = c(0, 1), 
-        expand = expansion(mult = c(0, 0)), 
-        breaks = pretty_breaks()) +  
-      # Prettify y axis...
-      scale_y_continuous(
-        name   = "Density", 
-        limits = c(0, 1), 
-        expand = expansion(mult = c(0, 0.1)), 
-        breaks = pretty_breaks())
-    
-    # Prettify theme
-    g = g + theme_classic() + 
-      theme(axis.text     = element_text(size = 10),
-            axis.title.x  = element_text(
-              size = 20, margin = margin(b = 10, t = 20)),
-            axis.title.y  = element_text(
-              size = 20, margin = margin(l = 10, r = 20)),
-            axis.line     = element_blank(),
-            strip.text    = element_text(size = 14),
-            strip.background = element_blank(), 
-            panel.border  = element_rect(
-              linewidth = 0.5, fill = NA),
-            panel.spacing = unit(1, "lines"),
-            legend.title  = element_blank(),
-            legend.text   = element_text(size = 14),
-            legend.key    = element_blank(),
-            legend.position = "right", 
-            legend.key.height = unit(2, "lines"),
-            legend.key.width  = unit(2, "lines"))
-    
-    # Save figure to file
-    if (d_v == "disease") save_fig(g, "X01")
-    if (d_v == "vaccine") save_fig(g, "X02")
-  }
-}
-
 # ---------------------------------------------------------
 # Plot age targets as defined by WIISE and VIMC coverage data
 # ---------------------------------------------------------
@@ -653,222 +473,6 @@ plot_coverage_age_density = function() {
   save_fig(g, "S08")
 }
 
-# **** External models ****
-
-# ---------------------------------------------------------
-# Plot outcomes from each external model
-# ---------------------------------------------------------
-plot_external_models = function() {
-  
-  extern_deaths_dt = table("extern_all_models") %>%
-    format_d_v_a_name()
-  
-  # ---- By scenario ----
-  
-  scenario_dt = extern_deaths_dt %>%
-    group_by(d_v_a_name, model, scenario, year) %>%
-    summarise(deaths = sum(deaths)) %>%
-    ungroup() %>%
-    group_by(d_v_a_name, model, scenario) %>%
-    mutate(cum_deaths = cumsum(deaths)) %>%
-    ungroup() %>%
-    as.data.table()
-  
-  g1 = ggplot(scenario_dt) +
-    aes(x = year, 
-        y = deaths, 
-        colour   = model, 
-        linetype = scenario) +
-    geom_line() +
-    facet_wrap(
-      facets = vars(d_v_a_name), 
-      scales = "free_y") +
-    scale_y_continuous(
-      labels = comma)
-  
-  g2 = ggplot(scenario_dt) +
-    aes(x = year, 
-        y = cum_deaths, 
-        colour   = model, 
-        linetype = scenario) +
-    geom_line() +
-    facet_wrap(
-      facets = vars(d_v_a_name), 
-      scales = "free_y") +
-    scale_y_continuous(
-      labels = comma)
-  
-  # ---- By impact ----
-  
-  impact_dt = scenario_dt %>%
-    group_by(d_v_a_name, model, year) %>%
-    mutate(averted = deaths[scenario == "no_vaccine"] - deaths) %>%
-    ungroup() %>%
-    filter(scenario == "vaccine") %>%
-    select(d_v_a_name, model, year, averted) %>%
-    group_by(d_v_a_name, model) %>%
-    mutate(cum_averted = cumsum(averted)) %>%
-    ungroup() %>%
-    as.data.table()
-  
-  g3 = ggplot(impact_dt) +
-    aes(x = year, 
-        y = averted, 
-        colour = model) +
-    geom_line() +
-    facet_wrap(
-      facets = vars(d_v_a_name), 
-      scales = "free_y") +
-    scale_y_continuous(
-      labels = comma)
-  
-  g4 = ggplot(impact_dt) +
-    aes(x = year, 
-        y = cum_averted, 
-        colour = model) +
-    geom_line() +
-    facet_wrap(
-      facets = vars(d_v_a_name), 
-      scales = "free_y") +
-    scale_y_continuous(
-      labels = comma)
-  
-  # ---- By region ----
-  
-  region_dt = extern_deaths_dt %>%
-    group_by(d_v_a_name, model, region, year) %>%
-    mutate(averted = deaths[scenario == "no_vaccine"] - deaths) %>%
-    ungroup() %>%
-    filter(scenario == "vaccine") %>%
-    select(d_v_a_name, model, region, year, averted) %>%
-    group_by(d_v_a_name, model, region) %>%
-    mutate(cum_averted = cumsum(averted)) %>%
-    ungroup() %>%
-    as.data.table()
-  
-  g5 = ggplot(region_dt) +
-    aes(x = year, 
-        y = averted, 
-        colour   = region, 
-        linetype = model) +
-    geom_line() +
-    facet_grid(
-      rows   = vars(d_v_a_name), 
-      cols   = vars(region),
-      scales = "free_y") +
-    scale_y_continuous(
-      labels = comma)
-  
-  g6 = ggplot(region_dt) +
-    aes(x = year, 
-        y = cum_averted, 
-        colour   = region, 
-        linetype = model) +
-    geom_line() +
-    facet_grid(
-      rows   = vars(d_v_a_name), 
-      cols   = vars(region),
-      scales = "free_y") +
-    scale_y_continuous(
-      labels = comma)
-  
-  save_fig(g1, "X05")
-  save_fig(g3, "X06")
-  save_fig(g5, "X07")
-}
-
-# ---------------------------------------------------------
-# Plot countries with missing coverage data
-# ---------------------------------------------------------
-plot_missing_data = function() {
-  
-  message("  - Plotting countries with missing coverage data")
-  
-  # Country population (most recent year)
-  pop_dt = table("wpp_pop") %>%
-    filter(year == max(o$years)) %>%
-    group_by(country) %>%
-    summarise(pop = sum(pop)) %>%
-    ungroup() %>%
-    mutate(pop_norm = sqrt(pop / max(pop))) %>%
-    as.data.table()
-  
-  # Identify countries with no data
-  plot_dt = table("coverage") %>%
-    left_join(y  = table("d_v_a"), 
-              by = "d_v_a_id") %>%
-    select(disease, country) %>%
-    unique() %>%
-    mutate(value = 0) %>%
-    # Wrangle to sparse, long format...
-    pivot_wider(id_cols     = "country", 
-                names_from  = "disease") %>%
-    pivot_longer(cols = -country, 
-                 names_to = "disease") %>%
-    # Append pop size for missing countries...
-    left_join(y  = pop_dt, 
-              by = "country") %>%
-    mutate(value = ifelse(is.na(value), pop_norm, NA)) %>%
-    # Append region details...
-    left_join(y  = table("country"), 
-              by = "country") %>%
-    select(disease, region, country = country_name, value) %>%
-    arrange(disease, region, country) %>%
-    # Set plotting order...
-    mutate(region  = fct_inorder(region),
-           country = fct_inorder(country)) %>% 
-    as.data.table()
-  
-  # Extract population limits for the data
-  limits  = c(0, max(plot_dt$value, na.rm = TRUE))
-  colours = colour_scheme("pals::brewer.reds", n = 10)
-  
-  # Plot missing data
-  g = ggplot(plot_dt) +
-    aes(x = disease, 
-        y = fct_rev(country), 
-        fill = value) + 
-    # Plot missing countries in colour...
-    geom_tile() + 
-    # Facet by region for readability...
-    facet_wrap(
-      facets = vars(region), 
-      scales = "free") + 
-    # Prettify y axis...
-    scale_y_discrete(
-      name = "Countries with missing data") +
-    # Set continuous colour bar...
-    scale_fill_gradientn(
-      na.value = "white",
-      colours  = colours,
-      limits   = limits, 
-      breaks   = limits, 
-      labels   = c("Smaller population", 
-                   "Larger population"))
-  
-  # Prettify theme
-  g = g + theme_classic() +
-    theme(axis.title.x  = element_blank(),
-          axis.title.y  = element_text(size = 28),
-          axis.text.x   = element_text(
-            size = 8, hjust = 1, angle = 50),
-          axis.text.y   = element_text(size = 8),
-          axis.line     = element_blank(),
-          panel.border  = element_rect(
-            linewidth = 0.5, fill = NA),
-          panel.spacing.x = unit(0.5, "lines"), 
-          panel.spacing.y = unit(0.5, "lines"), 
-          strip.text    = element_text(size = 12),
-          strip.background = element_blank(), 
-          legend.title  = element_blank(),
-          legend.text   = element_text(size = 10),
-          legend.position = "bottom", 
-          legend.key.width = unit(4, "cm"))
-  
-  # Save figure to file
-  save_fig(g, "X08")
-}
-
 # **** Static models ****
 
 # ---------------------------------------------------------
@@ -970,88 +574,6 @@ plot_gbd_estimates = function() {
   save_fig(g, "S12")
 }
 
-# ---------------------------------------------------------
-# Plot proportion of GBD burden we have coverage data for
-# ---------------------------------------------------------
-plot_gbd_missing = function() {
-  
-  message("  - Plotting GBD burden by vaccine coverage status")
-  
-  # Dictionary for status groups
-  status_dict = c(
-    trivial     = "Countries with NO coverage data",
-    non_trivial = "Countries with coverage data")
-  
-  # Countries which have non-trivial status
-  coverage_dt = table("coverage") %>%
-    left_join(y  = table("d_v_a"), 
-              by = "d_v_a_id") %>%
-    filter(source == "static") %>%
-    select(disease, country) %>%
-    unique() %>%
-    mutate(status = "non_trivial")
-  
-  # GBD disease burden by status group
-  status_dt = table("gbd_estimates") %>%
-    group_by(disease, country, year) %>%
-    summarise(deaths = sum(deaths_disease)) %>%
-    ungroup() %>%
-    left_join(y  = coverage_dt, 
-              by = c("disease", "country")) %>%
-    replace_na(list(status = "trivial")) %>%
-    left_join(y  = table("disease_name"), 
-              by = "disease") %>%
-    group_by(disease_name, year, status) %>%
-    summarise(deaths = sum(deaths)) %>%
-    ungroup() %>%
-    mutate(status = recode(status, !!!status_dict), 
-           status = factor(status, status_dict)) %>%
-    as.data.table()
-  
-  # Plot burden over time by vaccine coverage status
-  g = ggplot(status_dt) +
-    aes(x = year, 
-        y = deaths, 
-        fill = status) +
-    geom_bar(stat = "identity") +
-    # Facet by disease...
-    facet_wrap(
-      facets = vars(disease_name), 
-      scales = "free_y") + 
-    # Set colour scheme...
-    scale_fill_manual(
-      name   = "Age group",
-      values = c("red", "gray")) +
-    # Prettify y axis...
-    scale_y_continuous(
-      name   = "GBD-estimated number of deaths", 
-      labels = comma,
-      expand = expansion(mult = c(0, 0.05)), 
-      breaks = pretty_breaks())
-  
-  # Prettify theme
-  g = g + theme_classic() + 
-    theme(axis.title.x  = element_blank(),
-          axis.title.y  = element_text(
-            size = 20, margin = margin(l = 10, r = 20)),
-          axis.text     = element_text(size = 10),
-          axis.line     = element_blank(),
-          strip.text    = element_text(size = 14),
-          strip.background = element_blank(), 
-          panel.border  = element_rect(
-            linewidth = 0.5, fill = NA),
-          panel.spacing = unit(1, "lines"),
-          legend.title  = element_blank(),
-          legend.text   = element_text(size = 14),
-          legend.key    = element_blank(),
-          legend.position = "bottom", 
-          legend.key.height = unit(2, "lines"),
-          legend.key.width  = unit(2, "lines"))
-  
-  # Save to file
-  save_fig(g, "X09")
-}
-
 # ---------------------------------------------------------
 # Plot static model pathogen vaccine efficacy with immunity decay
 # ---------------------------------------------------------
@@ -1618,230 +1140,8 @@ plot_impute_perform = function(metric) {
   save_fig(g, paste0("S16", save_sub))
 }
 
-#-------------------------------------------------
-# Plot model choice per region
-#-------------------------------------------------
-plot_model_choice = function(metric) {
-  
-  message("  - Plotting model choice by region")
-  
-  # ---- Load results from imputation ----
-  
-  # Function to load imputation results
-  load_results_fn = function(id) {
-    
-    # Load file and extract model details
-    result = try_load(
-      pth  = o$pth$impute, 
-      file = paste1("impute", metric, id), 
-      throw_error = FALSE) %>%
-      pluck("model")
-    
-    return(result)
-  }
-  
-  # Load imputation results for all d-v-a
-  choice_dt = table("d_v_a") %>%
-    filter(source == "vimc") %>%
-    pull(d_v_a_id) %>%
-    lapply(load_results_fn) %>%
-    rbindlist() 
-  
-  # Group model choice by income 
-  income_choice_dt = choice_dt %>%
-    group_by(d_v_a_id) %>%
-    count(income, model_id) %>%
-    group_by(model_id) %>%         
-    mutate(prop = prop.table(n)) 
-  
-  # Group model choice by region 
-  region_choice_dt = choice_dt %>%
-    group_by(d_v_a_id) %>%
-    count(region, model_id) %>%
-    group_by(model_id) %>%         
-    mutate(prop = prop.table(n)) 
-  
-  # Plot model choice by region
-  g1 = ggplot(region_choice_dt) + 
-    aes(x = model_id,
-        y = n,
-        fill = region) + 
-    geom_bar(stat = "identity") + 
-    # Simple faceting with wrap labelling...
-    facet_wrap(
-      facets   = vars(d_v_a_id), 
-      labeller = label_wrap_gen(width = 30), 
-      scales   = "free") 
-  
-  # Prettify theme
-  g1 = g1 + theme_classic() + 
-    theme(axis.text     = element_text(size = 10),
-          axis.title.x  = element_text(
-            size = 20, margin = margin(b = 10, t = 20)),
-          axis.title.y  = element_text(
-            size = 20, margin = margin(l = 10, r = 20)),
-          axis.line     = element_blank(),
-          strip.text    = element_text(size = 14),
-          strip.background = element_blank(), 
-          panel.border  = element_rect(
-            linewidth = 1, fill = NA),
-          panel.spacing = unit(1, "lines"),
-          legend.title  = element_blank(),
-          legend.text   = element_text(size = 14),
-          legend.key    = element_blank(),
-          legend.position = "right", 
-          legend.key.height = unit(2, "lines"),
-          legend.key.width  = unit(2, "lines"))
-  
-  
-  # Plot model choice by income level
-  g2 = ggplot(income_choice_dt) + 
-    aes(x = model_id,
-        y = n,
-        fill = income) + 
-    geom_bar(stat = "identity") + 
-    # Simple faceting with wrap labelling...
-    facet_wrap(
-      facets   = vars(d_v_a_id), 
-      labeller = label_wrap_gen(width = 30), 
-      scales   = "free") 
-  
-  # Prettify theme
-  g2 = g2 + theme_classic() + 
-    theme(axis.text     = element_text(size = 10),
-          axis.title.x  = element_text(
-            size = 20, margin = margin(b = 10, t = 20)),
-          axis.title.y  = element_text(
-            size = 20, margin = margin(l = 10, r = 20)),
-          axis.line     = element_blank(),
-          strip.text    = element_text(size = 14),
-          strip.background = element_blank(), 
-          panel.border  = element_rect(
-            linewidth = 1, fill = NA),
-          panel.spacing = unit(1, "lines"),
-          legend.title  = element_blank(),
-          legend.text   = element_text(size = 14),
-          legend.key    = element_blank(),
-          legend.position = "right", 
-          legend.key.height = unit(2, "lines"),
-          legend.key.width  = unit(2, "lines"))
-  
-  # Save to file
-  save_fig(g1, "X10")
-  save_fig(g2, "X11")
-}
-
 # **** Impact functions ****
 
-# ---------------------------------------------------------
-# Exploratory plots of data used to fit impact functions
-# ---------------------------------------------------------
-plot_impact_data = function(metric) {
-  
-  message("  - Plotting impact function fitting data")
-  
-  # Load data used for impact function fitting
-  data_dt = read_rds("impact", "impact", metric, "data") %>%
-    format_d_v_a_name()
-  
-  # ---- Plot 1: impact per FVP over time ----
-  
-  # Impact per FVP over time
-  g1 = ggplot(data_dt) +
-    aes(x = year, 
-        y = impact_fvp, 
-        colour = country) +
-    geom_line(show.legend = FALSE) +
-    # Faceting with wrap labelling...
-    facet_wrap(
-      facets   = vars(d_v_a_name), 
-      labeller = label_wrap_gen(width = 30), 
-      scales   = "free_y") + 
-    # Set colour scheme...
-    scale_colour_manual(
-      values = colour_scheme(
-        map = "viridis::viridis", 
-        n   = n_unique(all_countries()))) +
-    # Prettify x axis...
-    scale_x_continuous(
-      expand = c(0, 0), 
-      breaks = pretty_breaks()) +  
-    # Prettify y axis...
-    scale_y_continuous(
-      name   = "Impact per fully vaccinated", 
-      labels = scientific,
-      limits = c(0, NA),
-      expand = expansion(mult = c(0, 0.05)), 
-      breaks = pretty_breaks())
-  
-  # Prettify theme
-  g1 = g1 + theme_classic() + 
-    theme(axis.text     = element_text(size = 8),
-          axis.text.x   = element_text(hjust = 1, angle = 50),
-          axis.title.x  = element_text(
-            size = 16, margin = margin(l = 10, r = 20)),
-          axis.title.y  = element_text(
-            size = 16, margin = margin(b = 10, t = 20)),
-          axis.line     = element_blank(),
-          strip.text    = element_text(size = 10),
-          strip.background = element_blank(), 
-          panel.border  = element_rect(
-            linewidth = 0.5, fill = NA),
-          panel.spacing = unit(0.5, "lines"))
-  
-  # Save figure to file
-  save_fig(g1, "X12")
-  
-  # ---- Plot 2: cumulative impact vs cumulative FVP ----
-  
-  # Cumulative FVPs vs cumulative deaths averted
-  g2 = ggplot(data_dt) + 
-    aes(x = fvps, 
-        y = impact, 
-        colour = country) +
-    geom_line(show.legend = FALSE) +
-    # Faceting with wrap labelling...
-    facet_wrap(
-      facets   = vars(d_v_a_name), 
-      labeller = label_wrap_gen(width = 30), 
-      scales   = "free") + 
-    # Set colour scheme...
-    scale_colour_manual(
-      values = colour_scheme(
-        map = "viridis::viridis", 
-        n   = n_unique(all_countries()))) +
-    # Prettify x axis...
-    scale_x_continuous(
-      name   = "Cumulative FVPs per capita (including new birth cohorts)", 
-      limits = c(0, NA),
-      expand = expansion(mult = c(0, 0.05)), 
-      breaks = pretty_breaks()) +  
-    # Prettify y axis...
-    scale_y_continuous(
-      name   = "Cumulative impact per capita (including new birth cohorts)", 
-      labels = scientific,
-      limits = c(0, NA),
-      expand = expansion(mult = c(0, 0.05)), 
-      breaks = pretty_breaks())
-  
-  # Prettify theme
-  g2 = g2 + theme_classic() + 
-    theme(axis.text     = element_text(size = 8),
-          axis.title.x  = element_text(
-            size = 16, margin = margin(l = 10, r = 20)),
-          axis.title.y  = element_text(
-            size = 16, margin = margin(b = 10, t = 20)),
-          axis.line     = element_blank(),
-          strip.text    = element_text(size = 10),
-          strip.background = element_blank(), 
-          panel.border  = element_rect(
-            linewidth = 0.5, fill = NA),
-          panel.spacing = unit(0.5, "lines"))
-  
-  # Save figure to file
-  save_fig(g2, "X13")
-}
-
 # ---------------------------------------------------------
 # Plot impact function evaluation
 # ---------------------------------------------------------
@@ -2083,217 +1383,34 @@ plot_model_selection = function(metric) {
               axis.title.y  = element_text(
                 size = 20, margin = margin(l = 10, r = 20)),
               axis.line     = element_blank(),
-              panel.border  = element_rect(
-                linewidth = 0.5, fill = NA),
-              legend.title  = element_blank(),
-              legend.text   = element_text(size = 12),
-              legend.key    = element_blank(),
-              legend.position = "right", 
-              legend.key.height = unit(2, "lines"),
-              legend.key.width  = unit(2, "lines"))
-    }
-    
-    return(g)
-  }
-  
-  # ---- A variety of plots ----
-  
-  # Create a variety of plots
-  g = list(
-    
-    # Plot by disease-vaccine-activity
-    a = plot_selection("d_v_a_name", stat = "proportion"), 
-    b = plot_selection("d_v_a_name", stat = "number"),
-    
-    # Plot by country
-    c = plot_selection("country", type = "density"))
-  
-  # Save figures of interest
-  save_sub = letters[which(metric %in% o$metrics)]
-  save_fig(g$a, paste0("S18", save_sub))
-}
-
-# ---------------------------------------------------------
-# Plot impact ratios - either all or initial only
-# ---------------------------------------------------------
-plot_impact_fvps = function(metric, scope) {
-  
-  # Set distance between vaccines
-  width = 0.3
-  
-  # Function for averaging (mean or median)
-  avg_fn = get("mean")
-  
-  # ---- Load data based on scope ----
-  
-  # All time plot
-  if (scope == "all_time") {
-    
-    message("  - Plotting all-time impact per FVP: ", metric)
-    
-    # Load initial ratio data
-    impact_dt = read_rds("history", "impact", metric, "data")
-    
-    # Set a descriptive y-axis title
-    y_lab = "Impact per fully vaccinated person (log10 scale)"
-  }
-  
-  # Initial year plot
-  if (scope == "initial") {
-    
-    message("  - Plotting initial impact per FVP: ", metric)
-    
-    # Load initial ratio data
-    impact_dt = read_rds("history", "initial_ratio", metric) %>%
-      rename(impact_fvp = initial_ratio)
-    
-    # Set a descriptive y-axis title
-    y_lab = "Initial impact per FVP used for back projection (log10 scale)"
-  }
-  
-  # ---- Classify by income status ----
-  
-  # Load income status dictionary
-  income_dict = table("income_dict")
-  
-  # Load income status of each country
-  income_dt = table("income_status") %>%
-    filter(year == max(year)) %>%
-    left_join(y  = income_dict, 
-              by = "income") %>%
-    select(country, income = income_name)
-  
-  # Classify by income group
-  data_dt = impact_dt %>%
-    filter(impact_fvp > 1e-12) %>%
-    # Append d_v_a description...
-    format_d_v_a_name() %>%
-    # Append income status description...
-    left_join(y  = income_dt, 
-              by = "country") %>%
-    mutate(income = factor(
-      x      = income, 
-      levels = rev(income_dict$income_name))) %>%
-    select(d_v_a = d_v_a_name, income, impact_fvp) %>%
-    unique()
-  
-  # ---- Plotting coordinates ----
-  
-  # Set x values for each d_v_a
-  x_major = data_dt %>%
-    group_by(d_v_a) %>%
-    summarise(order = avg_fn(impact_fvp)) %>%
-    ungroup() %>%
-    arrange(desc(order)) %>%
-    mutate(x_major = 1 : n()) %>%
-    select(-order) %>%
-    as.data.table()
-  
-  # Offset each income group
-  x_minor = data_dt %>%
-    select(income) %>%
-    unique() %>%
-    arrange(income) %>%
-    mutate(x_minor = seq(
-      from = -width / 2, 
-      to   = width / 2, 
-      length.out = n()))
-  
-  # Width of offset of points within income group
-  x_spray = width / (nrow(income_dict) * 2)
-  
-  # Bring all x and y values together
-  plot_dt = data_dt %>%
-    left_join(y  = x_major, 
-              by = "d_v_a") %>%
-    left_join(y  = x_minor, 
-              by = "income") %>%
-    mutate(x = x_major + x_minor) %>%
-    select(d_v_a, income, x, y = impact_fvp) %>%
-    arrange(x)
-  
-  # Extract bounds (after transformation)
-  lb = floor(min(log10(data_dt$impact_fvp)))
-  ub = ceiling(max(log10(data_dt$impact_fvp)))
-  
-  # ---- Vaccine and income status averages ----
-  
-  # Average by vaccine and income group
-  income_average_dt = plot_dt %>%
-    group_by(income, x) %>%
-    summarise(y = avg_fn(y)) %>%
-    ungroup() %>%
-    arrange(x) %>%
-    as.data.table()
-  
-  # Average by vaccine (over all income groups)
-  vaccine_average_dt = data_dt %>%
-    group_by(d_v_a) %>%
-    summarise(y = avg_fn(impact_fvp)) %>%
-    ungroup() %>%
-    left_join(y  = x_major,
-              by = "d_v_a") %>%
-    select(d_v_a, x = x_major, y) %>%
-    arrange(x) %>%
-    as.data.table()
-  
-  # ---- Produce primary plot ----
-  
-  # Plot impact per FVP
-  g = ggplot(plot_dt) +
-    aes(x = x, y = y) +
-    # Plot all points by vaccine and income...
-    geom_point(
-      mapping  = aes(colour = income),
-      alpha    = 0.1,
-      shape    = 16,
-      stroke   = 0,
-      position = position_jitter(
-        width = x_spray,
-        seed  = 1)) + 
-    # Average of each vaccine...
-    geom_point(
-      data   = vaccine_average_dt,
-      colour = "black",
-      shape = 95, # Horizontal lines
-      size  = 20) +
-    # Average of each income group...
-    geom_point(
-      data    = income_average_dt, 
-      mapping = aes(fill = income), 
-      color   = "black", 
-      shape   = 23, 
-      size    = 3) + 
-    # Prettify x axis...
-    scale_x_continuous(
-      breaks = x_major$x_major, 
-      labels = x_major$d_v_a) +
-    # Prettify y axis...
-    scale_y_continuous(
-      name   = y_lab,
-      trans  = "log10",
-      labels = scientific, 
-      limits = c(10 ^ lb, 10 ^ ub),
-      expand = c(0, 0),
-      breaks = 10 ^ rev(ub : lb))
+              panel.border  = element_rect(
+                linewidth = 0.5, fill = NA),
+              legend.title  = element_blank(),
+              legend.text   = element_text(size = 12),
+              legend.key    = element_blank(),
+              legend.position = "right", 
+              legend.key.height = unit(2, "lines"),
+              legend.key.width  = unit(2, "lines"))
+    }
+    
+    return(g)
+  }
   
-  # ---- Prettify and save ----
+  # ---- A variety of plots ----
   
-  # Prettify theme
-  g = g + theme_classic() + 
-    theme(axis.text     = element_text(size = 10),
-          axis.text.x   = element_text(hjust = 1, angle = 50),
-          axis.title.x  = element_blank(),
-          axis.title.y  = element_text(
-            size = 16, margin = margin(l = 10, r = 20)),
-          axis.line     = element_blank(),
-          panel.grid.major.y = element_line(linewidth = 0.25),
-          panel.border  = element_rect(
-            linewidth = 0.5, fill = NA),
-          panel.spacing = unit(1, "lines"))
+  # Create a variety of plots
+  g = list(
+    
+    # Plot by disease-vaccine-activity
+    a = plot_selection("d_v_a_name", stat = "proportion"), 
+    b = plot_selection("d_v_a_name", stat = "number"),
+    
+    # Plot by country
+    c = plot_selection("country", type = "density"))
   
-  # Save figure to file
-  save_fig(g, "X14")
+  # Save figures of interest
+  save_sub = letters[which(metric %in% o$metrics)]
+  save_fig(g$a, paste0("S18", save_sub))
 }
 
 # **** Historical impact ****
@@ -2713,8 +1830,8 @@ plot_temporal_impact = function(metric) {
           legend.key.height = unit(2, "lines"),
           legend.key.width  = unit(2, "lines"))
   
-  if (metric == "deaths") save_fig(g, "S03")
-  if (metric == "dalys")  save_fig(g, "S04")
+  if (metric == "deaths") save_fig(g, "S04")
+  if (metric == "dalys")  save_fig(g, "S05")
 }
 
 # ---------------------------------------------------------
@@ -3158,7 +2275,7 @@ plot_mortality_change = function() {
           legend.position = "none")
   
   # Save figure to file
-  save_fig(g, "S05")
+  save_fig(g, "S03")
 }
 
 # ---------------------------------------------------------
@@ -3345,498 +2462,6 @@ plot_survival_increase = function(log_age = FALSE) {
   save_fig(g, "3")
 }
 
-# ---------------------------------------------------------
-# Plot absolute and relative probability of death in 2024
-# ---------------------------------------------------------
-plot_prob_death_age = function() {
-  
-  message("  - Plotting probability of death by age")
-  
-  # Number of log2 age bins to plot
-  age_bins = 4
-  
-  # Dictionary for each vaccine case
-  #
-  # NOTE: Using reverse order here for prettier plot
-  case_dict = list( 
-    no_vaccine = "No historical vaccination",
-    vaccine    = "Vaccination as obserevd")
-  
-  # Construct a somewhat elaborate y axis label
-  y_label = bquote(
-    ~.("Marginal probability of death in n")
-    ^.("th") 
-    ~.("year of life")
-    ~.(paste0("(", max(o$years), ")")))
-  
-  # Deaths averted by vaccination
-  averted_dt = read_rds("history", "burden_averted_deaths") %>%
-    filter(year == max(o$years)) %>%
-    # Apply age structure...
-    left_join(y  = table("impact_age_multiplier"), 
-              by = "d_v_a_id", 
-              relationship = "many-to-many") %>%
-    mutate(impact = impact * scaler) %>%
-    # Summarise over region...
-    append_region_name() %>%
-    group_by(region, age) %>%
-    summarise(averted = sum(impact)) %>%
-    ungroup() %>%
-    # Avoid zero for log scaling ...
-    mutate(age = age + 1) %>%
-    filter(age %in% 2 ^ (0 : age_bins)) %>%
-    select(region, age, averted) %>%
-    as.data.table()
-  
-  # Probability of death by age
-  plot_dt = table("wpp_pop") %>%
-    filter(year == max(o$years)) %>%
-    left_join(y  = table("wpp_deaths"), 
-              by = c("country", "year", "age")) %>%
-    # Avoid zero for log scaling ...
-    mutate(age = age + 1) %>%
-    filter(age %in% 2 ^ (0 : age_bins)) %>%
-    # Average prob of death by region and year...
-    append_region_name() %>%
-    group_by(region, age) %>%
-    summarise(pop    = sum(pop), 
-              deaths = sum(deaths)) %>%
-    ungroup() %>%
-    # Append deaths averted...
-    left_join(y  = averted_dt, 
-              by = c("region", "age")) %>%
-    replace_na(list(averted = 0)) %>%
-    mutate(vaccine    = deaths / pop, 
-           no_vaccine = (deaths + averted) / pop) %>%
-    # Melt to long format ready for plotting...
-    select(region, age, vaccine, no_vaccine) %>%
-    pivot_longer(cols = c(vaccine, no_vaccine), 
-                 names_to = "case") %>%
-    # Set regional order...
-    arrange(desc(value)) %>%
-    mutate(region = fct_inorder(region)) %>%
-    # Set vaccine case order...
-    mutate(case = recode(case, !!!case_dict), 
-           case = factor(case, case_dict)) %>%
-    arrange(region, case, age) %>%
-    as.data.table()
-  
-  # Plot both scenarios as side-by-side bars
-  g = ggplot(plot_dt) +
-    aes(x = age, 
-        y = value, 
-        fill = case) +
-    geom_col(position = "dodge") +
-    # Facet by region...
-    facet_wrap(
-      facets = vars(region), 
-      scales = "free_y") +
-    # Set colour scheme...
-    scale_fill_manual(
-      values = colours_who("logo", n = 2)) + 
-    # Prettify x axis...
-    scale_x_continuous(
-      name   = "Age (log2 scale)",
-      trans  = "log2", 
-      breaks = 2 ^ (0 : age_bins)) +
-    # Prettify y axis...
-    scale_y_continuous(
-      name   = y_label,
-      labels = percent, 
-      limits = c(0, NA), 
-      expand = expansion(mult = c(0, 0.05)), 
-      breaks = pretty_breaks())
-  
-  # Prettify theme
-  g = g + theme_classic() + 
-    theme(axis.title   = element_text(size = 20),
-          axis.title.x = element_text(
-            margin = margin(t = 20, b = 10)),
-          axis.title.y = element_text(
-            margin = margin(l = 10, r = 20)),
-          axis.text    = element_text(size = 12),
-          axis.ticks   = element_blank(), 
-          axis.line    = element_line(linewidth = 0.25),
-          strip.text   = element_text(size = 18),
-          strip.background = element_blank(), 
-          panel.spacing = unit(1, "lines"), 
-          panel.grid.major.y = element_line(linewidth = 0.25),
-          legend.title = element_blank(),
-          legend.text  = element_text(size = 16),
-          legend.key   = element_blank(),
-          legend.position = "bottom", 
-          legend.key.height = unit(2, "lines"),
-          legend.key.width  = unit(2, "lines"))
-  
-  # Save figure to file
-  save_fig(g, "X15")
-}
-
-# ---------------------------------------------------------
-# Plot measles deaths in context of all cause deaths
-# ---------------------------------------------------------
-plot_measles_in_context = function() {
-  
-  message("  - Plotting measles in context")
-  
-  # Upper age bound for estimates
-  age_bound = 5
-  
-  # Metrics to plot
-  metric_dict = list(
-    cum  = "Cumulative number of deaths in children under %i",
-    abs  = "Annual number of deaths in children under %i",
-    norm = "Normalised cause of death for children under %i",
-    cov  = "Measles vaccine coverage (%i-%i)")
-  
-  # Cause of death
-  cause_dict = list(
-    all_cause = "All cause deaths",
-    measles   = "Measles as cause of death", 
-    other     = "Other cause of death")
-  
-  # Define vaccine scenarios
-  case_dict = list( 
-    vaccine    = "Vaccination as obserevd", 
-    no_vaccine = "No historical vaccination")
-  
-  # Name of d-v-a to plot
-  measles_id = "Measles"
-  
-  # ID of coverage values to plot
-  coverage_ids = c(101, 102)
-  
-  # Colour schemes
-  colours = list(
-    area  = c("#B0B0B0", "#0189C5"), 
-    line  = c("#000000", "#0189C5"), 
-    doses = c("#C70A0A", "#901616", "#400808"))
-  
-  # ---- Deaths by cause ----
-  
-  # Death data - WPP and external measles models
-  deaths = list(
-    wpp    = table("wpp_deaths"), 
-    extern = table("extern_deaths"))
-  
-  # Measles deaths with and without vaccine
-  measles_deaths = deaths$extern %>%
-    format_d_v_a_name() %>%
-    # Pathogen and age group of interest...
-    filter(d_v_a_name == measles_id, 
-           age        <= age_bound) %>%
-    pivot_longer(cols = c(vaccine, no_vaccine), 
-                 names_to = "case") %>%
-    # Summarise over age bins of interest...
-    group_by(case, year) %>%
-    summarise(measles = sum(value)) %>%
-    ungroup() %>%
-    as.data.table()
-  
-  # All cause deaths from WPP
-  all_deaths = deaths$wpp %>%
-    # Age group of interest...
-    filter(age <= age_bound) %>%
-    group_by(year) %>%
-    summarise(deaths = sum(deaths)) %>%
-    ungroup() %>%
-    # Repeat for each vaccine scenario...
-    expand_grid(case = names(case_dict)) %>%
-    select(case, year, deaths) %>%
-    # Append number of deaths averted from vaccination...
-    left_join(y  = measles_deaths, 
-              by = c("case", "year")) %>%
-    group_by(year) %>%
-    mutate(averted = measles - measles[case == "vaccine"]) %>%
-    ungroup() %>%
-    # Increment all cause deaths in scenario of no vaccination...
-    mutate(all_cause = deaths + averted) %>%
-    select(case, year, all_cause) %>%
-    arrange(case, year) %>%
-    as.data.table()
-  
-  # Combine causes of death for measles in context
-  context_dt = all_deaths %>%
-    left_join(y  = measles_deaths, 
-              by = c("case", "year")) %>%
-    mutate(other = all_cause - measles) %>% 
-    # Melt cause of death to long format...
-    pivot_longer(cols = names(cause_dict), 
-                 names_to  = "cause", 
-                 values_to = "abs") %>%
-    # Calculate normalised cause of death...
-    group_by(case, year) %>%
-    mutate(norm = abs / abs[cause == "all_cause"]) %>%
-    ungroup() %>%
-    # Calculate cumulative values...
-    group_by(case, cause) %>%
-    mutate(cum = cumsum(abs)) %>%
-    ungroup() %>%
-    # Convert to tidy format...
-    pivot_longer(cols = any_names(metric_dict), 
-                 names_to = "metric") %>%
-    # Set appropriate order...
-    select(case, metric, cause, year, value) %>%
-    arrange(case, metric, cause, year) %>%
-    as.data.table()
-  
-  # ---- Construct plotting datatables ----
-  
-  # First update metric dictionary with values
-  metric_dict = list(
-    cum  = sprintf(metric_dict$cum,  age_bound),
-    abs  = sprintf(metric_dict$abs,  age_bound),
-    norm = sprintf(metric_dict$norm, age_bound),
-    cov  = sprintf(metric_dict$cov,  min(o$years), max(o$years)))
-  
-  # Lines for measles and all cause
-  line_dt = context_dt %>%
-    filter(cause != "other") %>%
-    # Don't plot normalised all cause as always one...
-    filter(!(metric == "norm" & 
-               cause == "all_cause")) %>%
-    # Recode variables and set ordering...
-    mutate(case   = recode(case,   !!!case_dict), 
-           case   = factor(case,   case_dict), 
-           metric = recode(metric, !!!metric_dict), 
-           metric = factor(metric, metric_dict), 
-           cause  = recode(cause,  !!!cause_dict),
-           cause  = factor(cause,  cause_dict))
-  
-  # Area for measles and other cause in vaccine case only
-  area_dt = context_dt %>%
-    filter(cause != "all_cause", 
-           case  == "vaccine") %>%
-    select(-case) %>%
-    # Recode variables and set ordering...
-    mutate(metric = recode(metric, !!!metric_dict), 
-           metric = factor(metric, metric_dict), 
-           cause  = recode(cause,  !!!cause_dict),
-           cause  = factor(cause,  cause_dict))
-  
-  # Gloabl measles coverage over time
-  coverage_dt = table("coverage_global") %>%
-    filter(d_v_a_id %in% coverage_ids) %>%
-    # Append metric details...
-    left_join(y  = table("vaccine_name"), 
-              by = "vaccine") %>%
-    mutate(metric = metric_dict$cov, 
-           metric = factor(metric, metric_dict)) %>%
-    # Tidy up...
-    select(metric, vaccine_name, year, value = coverage) %>%
-    as.data.table()
-  
-  # ---- Produce plot ----
-  
-  # Plot measles in context of all cause deaths
-  g = ggplot(area_dt) +
-    aes(x = year, 
-        y = value) +
-    # Plot cause of death...
-    geom_area(
-      mapping = aes(fill = fct_rev(cause)),
-      alpha   = 0.5,
-      show.legend = FALSE) +
-    # Plot all cause deaths...
-    geom_line(
-      data    = line_dt, # [case == case_dict[[1]], ]
-      mapping = aes(
-        colour = cause, 
-        linetype = case), 
-      linewidth = 1.1) +
-    # Set primary colour scheme...
-    scale_colour_manual(
-      name   = "Cause of death", 
-      values = colours$line) + 
-    scale_fill_manual(
-      values = colours$area) +
-    # Prettify primary legends...
-    guides(colour = guide_legend(
-      order   = 2, 
-      reverse = TRUE)) + 
-    # Plot number of doses...
-    ggnewscale::new_scale_color() + 
-    geom_line(
-      data    = coverage_dt,
-      mapping = aes(colour = vaccine_name)) +
-    # Set doses colour scheme...
-    scale_colour_manual(
-      name   = "Vaccine", 
-      values = colours$doses) + 
-    # Facet by metric...
-    facet_wrap(
-      facets = vars(metric),
-      scales = "free_y",
-      ncol   = 1) +
-    # Prettiy x axis...
-    scale_x_continuous(
-      limits = c(min(o$years), max(o$years)), 
-      expand = expansion(mult = c(0, 0)), 
-      breaks = seq(min(o$years), max(o$years), by = 5)) +
-    # Prettify y axis...
-    expand_limits(y = c(0, 1)) + 
-    scale_y_continuous(
-      labels = comma, 
-      expand = expansion(mult = c(0, 0)), 
-      breaks = pretty_breaks()) +
-    # Prettify primary legends...
-    guides(
-      colour    = guide_legend(
-        order   = 3), 
-      linetype  = guide_legend(
-        order   = 1, 
-        title   = "Vaccine scenario"))
-  
-  # Prettify theme
-  g = g + theme_classic() + 
-    theme(axis.title    = element_blank(),
-          axis.text     = element_text(size = 10),
-          axis.text.x   = element_text(hjust = 1, angle = 50), 
-          axis.line     = element_blank(),
-          strip.text    = element_text(size = 16),
-          strip.background = element_blank(), 
-          panel.border  = element_rect(
-            linewidth = 0.5, fill = NA),
-          panel.spacing = unit(1, "lines"),
-          legend.title  = element_text(size = 14),
-          legend.text   = element_text(size = 12),
-          legend.margin = margin(t = 40, b = 40, l = 10),
-          legend.position = "right", 
-          legend.key.height = unit(2, "lines"),
-          legend.key.width  = unit(2, "lines"))
-  
-  # Save figure to file
-  save_fig(g, "X16")
-}
-
-# ---------------------------------------------------------
-# Plot comparison of EPI50 outcomes vs VIMC outcomes
-# ---------------------------------------------------------
-plot_vimc_comparison = function() {
-  
-  message("  - Plotting comparison of EPI50 vs VIMC outcomes")
-  
-  # Dictionary for source of estimates
-  type_dict = c(
-    impact_epi50 = "Deaths averted: EPI50", 
-    impact_vimc  = "Deaths averted: VIMC", 
-    fvps_epi50   = "Fully vaccinated people: EPI50", 
-    fvps_vimc    = "Fully vaccinated people: VIMC")
-  
-  # Dictionary for metrics of interest
-  metric_dict = c(
-    impact = "Deaths averted", 
-    fvps   = "Fully vaccinated people")
-  
-  # ---- Load EPI50 and VIMC outcomes ----
-  
-  # Store all plotting values in list
-  plot_list = list(
-    
-    # EPI50 impact 
-    a = read_rds("history", "burden_averted_deaths") %>%
-      lazy_dt() %>%
-      group_by(d_v_a_id) %>%
-      summarise(value = round(sum(impact))) %>%
-      ungroup() %>%
-      mutate(metric = "impact", 
-             type   = "epi50") %>%
-      as.data.table(), 
-    
-    # VIMC impact
-    b = table("vimc_estimates") %>%
-      lazy_dt() %>%
-      group_by(d_v_a_id) %>%
-      summarise(value = sum(deaths_averted)) %>%
-      ungroup() %>%
-      mutate(metric = "impact", 
-             type   = "vimc") %>%
-      as.data.table(), 
-    
-    # EPI50 FVPs
-    c = table("coverage") %>%
-      lazy_dt() %>%
-      group_by(d_v_a_id) %>%
-      summarise(value = round(sum(fvps))) %>%
-      ungroup() %>%
-      mutate(metric = "fvps", 
-             type   = "epi50") %>%
-      as.data.table(),
-    
-    # VIMC FVPs
-    d = table("coverage_source") %>%
-      filter(source == "vimc") %>%
-      lazy_dt() %>%
-      group_by(d_v_a_id) %>%
-      summarise(value = round(sum(fvps))) %>%
-      ungroup() %>%
-      mutate(metric = "fvps", 
-             type   = "vimc") %>%
-      as.data.table())
-  
-  # ---- Concatenate all outcomes ---- 
-  
-  # Combine all plotting data
-  plot_dt = rbindlist(plot_list) %>%
-    left_join(y  = table("d_v_a"), 
-              by = "d_v_a_id") %>%
-    filter(source == "vimc") %>%
-    # Append d_v_a names...
-    select(d_v_a_name, metric, type, value) %>%
-    mutate(type = paste1(metric, type)) %>%
-    # Set plotting order...
-    mutate(metric = recode(metric, !!!metric_dict), 
-           metric = factor(metric, metric_dict), 
-           type   = recode(type,   !!!type_dict), 
-           type   = factor(type,   type_dict)) %>%
-    as.data.table()
-  
-  # ---- Produce plot ----
-  
-  # Plot bars of EPI50 outcomes vs VIMC
-  g = ggplot(plot_dt) + 
-    aes(x = d_v_a_name, 
-        y = value, 
-        fill = type) + 
-    geom_col(position = "dodge") + 
-    # Facet by metric...
-    facet_wrap(
-      facets = vars(metric), 
-      scales = "free_y", 
-      ncol   = 1) + 
-    # Set paired colour scheme...
-    scale_fill_manual(
-      values = rev(colour_scheme(
-        map = "brewer::paired", 
-        n   = length(type_dict)))) + 
-    # Prettify y axis...
-    scale_y_continuous(
-      labels = comma)
-  
-  # Prettify theme
-  g = g + theme_classic() + 
-    theme(axis.title    = element_blank(),
-          axis.text     = element_text(size = 11),
-          axis.text.x   = element_text(hjust = 1, angle = 50), 
-          axis.line     = element_blank(),
-          strip.text    = element_text(size = 14),
-          strip.background = element_blank(), 
-          panel.border  = element_rect(
-            linewidth = 0.5, fill = NA),
-          panel.spacing = unit(1, "lines"),
-          panel.grid.major.y = element_line(linewidth = 0.25),
-          legend.title  = element_blank(),
-          legend.text   = element_text(size = 11),
-          legend.key    = element_blank(),
-          legend.position = "right", 
-          legend.spacing.y  = unit(1, "lines"),
-          legend.key.height = unit(2, "lines"),
-          legend.key.width  = unit(2, "lines"))
-  
-  # Save theis figure to file
-  save_fig(g, "X17")
-}
-
 # **** Helper functions ****
 
 # ---------------------------------------------------------
@@ -3945,32 +2570,6 @@ format_d_v_a_name = function(dt) {
   return(order_dt)
 }
 
-# ---------------------------------------------------------
-# Get vector of colours for a given variable
-# ---------------------------------------------------------
-get_palette = function(variable) {
-  
-  # Input argument must be specified in o$palette
-  if (!variable %in% names(o$palette))
-    stop("Input '", variable, "' must be one of: ", 
-         paste(names(o$palette), collapse = ", "))
-  
-  # Number of colours to generate based on input argument
-  n_colours = list(
-    disease = n_unique(table("d_v_a")$disease), 
-    region  = n_unique(table("country")$region), 
-    income  = n_unique(table("income_status")$income))
-  
-  # Shorthand for palette and number of colours
-  p = o$palette[[variable]]
-  n = n_colours[[variable]]
-  
-  # Generate vector of colours
-  colours = colour_scheme(p, n = n)
-  
-  return(colours)
-}
-
 # ---------------------------------------------------------
 # Save a ggplot figure to file with default settings
 # ---------------------------------------------------------
diff --git a/regression.R b/regression.R
index 4b98770..bbabfeb 100644
--- a/regression.R
+++ b/regression.R
@@ -80,9 +80,6 @@ run_regression = function(case, metric) {
   
   # Plot predicted vs observed for each country
   plot_impute_perform(metric)
-  
-  # Plot model selection by disease
-  plot_model_choice(metric)
 }
 
 # ---------------------------------------------------------
diff --git a/results.R b/results.R
index 71949eb..2721c12 100644
--- a/results.R
+++ b/results.R
@@ -27,27 +27,11 @@ run_results = function() {
 
     # Total number of FVP over time by source
     plot_total_fvps()
-    plot_smooth_fvps()
-
-    # Plot coverage density by disease
-    plot_coverage()
     
     # Coverage data density by age
     plot_coverage_age_density()
   }
   
-  # ---- External model plots ----
-  
-  # Check plotting flag
-  if (o$plot_external) {
-    
-    # Plot outcomes from each external model
-    plot_external_models()
-    
-    # Missing coverage data by country
-    plot_missing_data()
-  }
-  
   # ---- Static model plots ----
   
   # Check plotting flag
@@ -55,10 +39,7 @@ run_results = function() {
     
     # Global Burden of Disease death estimates by age
     plot_gbd_estimates()
-
-    # Proportion of GBD burden we have coverage data for
-    plot_gbd_missing()
-
+    
     # Plot vaccine efficacy profiles for static model pathogens
     plot_vaccine_efficacy()
 
@@ -82,9 +63,6 @@ run_results = function() {
       
       # Plot predicted vs observed for each country
       plot_impute_perform(metric)
-      
-      # Plot model selection by disease
-      plot_model_choice(metric)
     }
   }
   
@@ -96,9 +74,6 @@ run_results = function() {
     # Repeat for deaths and DALYs
     for (metric in o$metrics) {
       
-      # Exploratory plots of data used to fit impact functions
-      plot_impact_data(metric)
-      
       # Plot impact function evaluation
       plot_model_fits(metric)
 
@@ -131,19 +106,6 @@ run_results = function() {
     
     # Plot absolute and relative probability of death in 2024
     plot_survival_increase()
-    
-    # Plot absolute and relative probability of death in 2024
-    plot_prob_death_age()
-    
-    # Measles deaths in context of all cause deaths
-    plot_measles_in_context()
-    
-    # Inital impact ratios used to back project
-    for (metric in o$metrics)
-      plot_impact_fvps(metric, scope = "initial")
-    
-    # Plot comparison of EPI50 outcomes vs VIMC outcomes
-    plot_vimc_comparison()
   }
   
   # ---- Main results table ----