Fixes to allow looking at BA.2 vs BA.1

R/estimate-generation-time.R

@@ -7,7 +7,7 @@ gt_load_model <- function(model = here::here("stan/generation-time.stan"),
 
 gt_prior <- function(type = "household", source = "hart2021") {
   type <- match.arg(type, choices = c("household", "intrinsic"))
-  source <- match.arg(source, choices = c("hart2021"))
+  source <- match.arg(source, choices = c("hart2021", "abbott2022"))
 
   if (type %in% "intrinsic" & source %in% "hart2021") {
     gt <- list(
@@ -21,7 +21,7 @@ gt_prior <- function(type = "household", source = "hart2021") {
       source = "hart2021",
       doi = "10.1101/2021.10.21.21265216v1"
     )
-  }else if (type %in% "household" & source %in% "hart2021") {
+  } else if (type %in% "household" & source %in% "hart2021") {
     gt <- list(
       # From Hart et al.
       # https://www.medrxiv.org/content/10.1101/2021.10.21.21265216v1
@@ -33,7 +33,32 @@ gt_prior <- function(type = "household", source = "hart2021") {
       source = "hart2021",
       doi = "10.1101/2021.10.21.21265216v1"
     )
+  } else if (type %in% "intrinsic" & source %in% "abbott2022") {
+    gt <- list(
+      # From Abbott et al.
+      # https://www.medrxiv.org/content/10.1101/2022.01.08.22268920v1
+      # Assuming symmetric normal which is incorrect but an approximation
+      mean_mean = 3.3,
+      mean_sd = 0.7,
+      sd_mean = 3.5,
+      sd_sd = 1.2,
+      source = "abbott2022",
+      doi = "10.1101/2022.01.08.22268920"
+    )
+  } else if (type %in% "household" & source %in% "abbott2022") {
+    gt <- list(
+      # From Abbott et al.
+      # https://www.medrxiv.org/content/10.1101/2022.01.08.22268920v1
+      # Assuming symmetric normal which is incorrect but an approximation
+      mean_mean = 2.2,
+      mean_sd = 0.3,
+      sd_mean = 2.7,
+      sd_sd = 1.1,
+      source = "abbott2022",
+      doi = "10.1101/2022.01.08.22268920"
+    )
   }
+
   return(gt)
 }
 

R/load-data.R

@@ -11,15 +11,14 @@ source(here::here("R", "load-public-data.R"))
 source(here::here("R", "load-private-data.R"))
 source(here::here("R", "munge-data.R"))
 
-load_data <- function(open_data = TRUE, data_type = "specimen",
-                      start_date = as.Date("2021-11-20")) {
+load_data <- function(open_data = TRUE, ...) {
   if (open_data) {
     if (data_type == "onset") {
       stop("Public data are by specimen date only")
     }
     regional <- load_public_data()
   } else {
-    regional <- load_private_data(data_type = data_type)
+    regional <- load_private_data(...)
   }
   regional <- add_england_totals(regional)
   return(regional)

R/load-private-data.R

@@ -7,25 +7,37 @@ library(zoo)
 library(ggplot2)
 
 
-load_private_date <- function(path = "data/private", data_type = "specimen") {
+load_private_data <- function(path = "data/private", id_type = "sgtf", data_type = "specimen", clean = TRUE, start_date = NULL) {
   data_path <- here(path,
-                    paste0("sgtf_",
-                           ifelse(data_type == "onset", data_type,
-                                  "daily"),
-                           "_england.csv"))
+                    paste0(id_type, "_",
+                           data_type,
+                           "_england",
+			   ifelse(clean, "_clean", ""),
+			   ".csv"))
 
   regional_raw <- readr::read_csv(data_path, show_col_types = FALSE)
 
+  if ("prop" %in% colnames(regional_raw)) {
+    regional_raw <- regional_raw %>%
+    select(-prop)
+  }
+  if (id_type == "omicron") {
+    regional_raw <- regional_raw %>%
+      rename(sgtf = Omicron, non_sgtf = `Omicron BA.2`)
+  }
   regional <- regional_raw %>%
-    select(-prop) %>%
     rename(date = grep("date", names(.), value = TRUE),
            region = nhser_name,
            sgtf_unknown = "NA") %>%
-    rowwise() %>%
-    mutate(total_sgt = sum(non_sgtf, sgtf, na.rm = TRUE),
-           total_cases = sum(non_sgtf, sgtf, sgtf_unknown, na.rm = TRUE)) %>%
-    ungroup() %>%
+    replace_na(list(sgtf_unknown = 0, non_sgtf = 0, sgtf = 0)) %>%
+    mutate(total_sgt = non_sgtf + sgtf,
+           total_cases = non_sgtf+ sgtf + sgtf_unknown) %>%
     mutate(region = tidyr::replace_na(region, "Unknown"),
            source = "private")
+
+  if (!is.null(start_date)) {
+    regional <- regional %>%
+      filter(date >= start_date)
+  }
   return(regional)
-}
+}

data/retrospective/posterior_summary.csv

@@ -1,17 +1,17 @@
 stratification,gt_type,variable,mean,median,sd,mad,q5,q95,rhat,ess_bulk,ess_tail
-region,household,gt_mean,3.225741125,3.222055,0.464544029530439,0.46309011,2.473732,3.988507,1.00080357267048,3251.09444256016,2547.51148832665
-region,household,gt_sd,2.423210335,2.426655,0.329790304735315,0.332361855,1.873326,2.978121,1.00164478915307,3354.5703472291,2449.16892989282
-region,household,voc_gt_mean_mod,0.68792100425,0.6514595,0.190666490343531,0.1499027208,0.46310245,1.0258795,1.00300278280105,1195.42661611848,1460.3464740799
-region,household,voc_gt_sd_mod,1.120743944525,1.079765,0.382727999766832,0.344682261,0.5654927,1.7711895,1.00434808199062,979.401527971105,1476.63998679927
-region,household,voc_gt_mean,2.192822175,2.11112,0.572569308051742,0.494321079,1.453528,3.211919,1.0033447166854,1350.30387518086,1643.01626315082
-region,household,voc_gt_sd,2.73216737375,2.591325,1.07806517522467,0.910516551,1.2846985,4.605492,1.00399595377353,1038.35400630956,1512.36496529864
-region,household,ta,1.84961721,1.839775,0.155671823916874,0.151625502,1.616737,2.122095,1.00123934901765,2864.69056327974,2762.81279840174
-region,household,sigma,0.0224171332162575,0.02265575,0.0119884223412889,0.01349329086,0.0030806955,0.041832715,1.00710736931297,367.308514343612,839.154396390054
-region,intrinsic,gt_mean,4.59910677,4.599845,0.36558537041214,0.358967112,3.99231,5.193273,1.00042598698629,6200.89377665651,3109.04819224535
-region,intrinsic,gt_sd,3.10512604,3.106155,0.172802975907503,0.175828947,2.8225685,3.3864905,1.00034361882293,6498.82151907565,3061.62451342061
-region,intrinsic,voc_gt_mean_mod,0.720927093,0.699979,0.160257431385769,0.1340678115,0.51530905,0.988511,1.00171590192607,1219.52666007364,1447.91751714838
-region,intrinsic,voc_gt_sd_mod,1.13295678275,1.09178,0.371678684630054,0.3188672298,0.6347775,1.7318365,1.00215535731066,1143.42192084155,1464.74025112686
-region,intrinsic,voc_gt_mean,3.30269248,3.215225,0.718781057119181,0.616361298,2.3629625,4.47496,1.00161843507025,1221.56085943347,1367.07214986368
-region,intrinsic,voc_gt_sd,3.523606806,3.40466,1.19323532260682,1.020273429,1.9320165,5.4527565,1.00134635891911,1129.92192761031,1281.44647493737
-region,intrinsic,ta,2.446883145,2.43633,0.178791349856767,0.175191429,2.1801225,2.753062,1.0018317860861,2619.75030654476,2088.52829056239
-region,intrinsic,sigma,0.020420386529275,0.0203062,0.0117707642371426,0.01359151311,0.002148792,0.04017923,1.00920944199813,467.064211939111,947.607178049898
+region,household,gt_mean,2.2646718125,2.269825,0.306006620109153,0.309307425,1.7443205,2.762074,1.00201369399288,2237.72829683861,822.968327647243
+region,household,gt_sd,1.807085904012,1.845865,0.748470720916594,0.788202051,0.4949534,2.998071,1.0063217181471,817.779694452818,1198.3177686017
+region,household,voc_gt_mean_mod,1.2216117995,1.14463,0.260159910565222,0.183508815,0.9558883,1.8025525,1.00919914194233,589.935128799332,313.224554120584
+region,household,voc_gt_sd_mod,1.11917299005,1.074005,0.627038187841528,0.670757892,0.2320322,2.1478505,1.00337146553337,1393.29078031312,882.394154505056
+region,household,voc_gt_mean,2.75808659,2.63217,0.654234394187325,0.563135958,1.9163275,4.0432875,1.00464718938323,1198.3326637478,1034.98384732952
+region,household,voc_gt_sd,2.130635955115,1.71571,1.64982883969248,1.561748601,0.22652175,5.635851,1.00705248165809,639.074313623767,476.88603736515
+region,household,ta,1.33568829,1.33118,0.062098204535055,0.0595856940000001,1.241469,1.447237,1.00086647586091,1854.60600465136,1517.91166744686
+region,household,sigma,0.00227917784682075,0.001885205,0.00177342146645277,0.0016987638213,0.00013828995,0.0058055615,1.00349872370038,1276.42421252123,746.614020544549
+region,intrinsic,gt_mean,3.597418045,3.584055,0.642684440154834,0.658593159,2.5640835,4.6573225,1.00373318264472,2251.23545963828,2638.69197971454
+region,intrinsic,gt_sd,2.447875084325,2.461765,0.86317987454556,0.871005261,0.9407665,3.838004,1.00394524460998,968.63394258405,458.599276872034
+region,intrinsic,voc_gt_mean_mod,1.24739859125,1.172805,0.259429139589357,0.191522268,0.97400755,1.7901885,1.00294087144331,1021.33690923067,479.740320843753
+region,intrinsic,voc_gt_sd_mod,1.100422840025,1.059765,0.585027499263111,0.6351562182,0.26189515,2.0816735,1.00304598607234,1339.25199203725,633.084884091115
+region,intrinsic,voc_gt_mean,4.4872535825,4.297415,1.25738964463698,1.032126816,2.8787545,6.8100775,1.00339369860237,971.833685978459,500.699663616734
+region,intrinsic,voc_gt_sd,2.8670008490275,2.398075,2.14161569765119,1.982888544,0.3819898,6.914455,1.00324239959122,979.270471987653,455.589016468813
+region,intrinsic,ta,1.597935975,1.582015,0.151905151691366,0.145976796,1.3770355,1.8698365,1.00254706519252,1143.08189843651,545.902799670362
+region,intrinsic,sigma,0.002273865474545,0.001919925,0.00171733380517185,0.0017201481024,0.00019400885,0.005587105,1.00021703703314,2796.39691384006,1910.68143849923

scripts/estimate-observed-transmission-parameters.R

@@ -16,7 +16,7 @@ source(here("R", "estimate-generation-time.R"))
 # Load growth estimates
 region_target_date <- get_latest_date()
 region_growth <- load_growth(
-  region_target_date, min_date = "2021-12-01", max_date = "2021-12-23"
+  region_target_date, min_date = "2022-01-12", max_date = "2022-02-06"
 )
 region_growth <- region_growth[!(region %in% "England")]
 
@@ -33,7 +33,7 @@ grid <- CJ(
 )
 
 grid[, gt_prior := purrr::map(
-  gt_type, ~ gt_prior(source = "hart2021", type = .x))
+  gt_type, ~ gt_prior(source = "abbott2022", type = .x))
 ]
 
 grid <- merge(grid, growth, by = "stratification")
@@ -71,4 +71,4 @@ fwrite(
 fwrite(
   posterior_predictions,
   here::here("data", "retrospective", "posterior_predictions.csv")
-)
+)

scripts/update-regional-data.R

@@ -4,11 +4,13 @@ library(ggplot2)
 source(here::here("R", "load-data.R"))
 source(here::here("R", "plot-daily-cases.R"))
 
-daily <- load_data()
+start_date <- as.Date("2021-11-15")
+
+daily <- load_data(open_data = FALSE, data_type = "onset", start_date = start_date, clean = FALSE)
 
 date <- max(daily$date)
 download_date <- Sys.Date()
 readr::write_csv(
   daily,
   here::here("data", "public", paste0(download_date, "-cases-by-sgtf.csv"))
-)
+)

scripts/update-regional-estimates.R

@@ -33,13 +33,17 @@ target_date <- get_latest_date()
 results <- load_results(target_date)
 
 # Estimation start date
-start_date <- as.Date("2021-11-16")
-start_sgtf_date <- as.Date("2021-11-23")
+start_date <- as.Date("2022-01-05")
+start_sgtf_date <- as.Date("2022-01-12")
 
 # Load data for the target date
 daily_regional <- load_local_data(target_date) %>%
   filter(date >= start_date)
 
+# for BA.2: swap sgtf and non-sgtf
+daily_regional <- daily_regional %>%
+  rename(sgtf = non_sgtf, non_sgtf = sgtf)
+
 # check is the sgtf data is newer than the sgtf data in the results
 # if results are present
 if (!is.null(results$posterior)) {
@@ -91,7 +95,7 @@ sgtf_regional <- daily_regional %>%
   filter(!(is.na(cases) & is.na(seq_voc)))
 
 # Estimate models for SGTF data
-region_omicorn_forecasts <- build_models_by_region(
+region_omicron_forecasts <- build_models_by_region(
   sgtf_regional, sgtf_parameters,
   variant_relationships = c("scaled", "correlated"),
   cores_per_model = 2, chains = 2, samples_per_chain = 2000,
@@ -101,14 +105,14 @@ region_omicorn_forecasts <- build_models_by_region(
 omicron_results <- list(
   data = sgtf_regional,
   posterior = summary(
-    region_omicorn_forecasts, target = "posterior", type = "all"
+    region_omicron_forecasts, target = "posterior", type = "all"
   )[,
    loo := NULL
   ],
-  diagnostics = summary(region_omicorn_forecasts, target = "diagnostics")[,
+  diagnostics = summary(region_omicron_forecasts, target = "diagnostics")[,
    loo := NULL
   ],
-  loo = extract_loo(region_omicorn_forecasts)
+  loo = extract_loo(region_omicron_forecasts)
 )
 
 save_results(omicron_results, "sgtf", target_date)

writeup/generation-time.Rmd

@@ -44,9 +44,9 @@ source(here("R", "estimate-generation-time.R"))
 ```
 
 ```{r setup-dates}
-date_data_start <- as.Date("2021-11-23")
-date_gt_start <- as.Date("2021-12-01")
-date_forecast_start <- as.Date("2021-12-23")
+date_data_start <- as.Date("2022-01-01")
+date_gt_start <- as.Date("2022-01-12")
+date_forecast_start <- as.Date("2022-02-06")
 ```
 
 1. Centre for the Mathematical Modelling of Infectious Diseases, London School of Hygiene & Tropical Medicine, London WC1E 7HT, United Kingdom
@@ -94,7 +94,7 @@ A full description of this model can be found in the documentation for the [`for
 
 ```{r load-regional-growth}
 # get regional estimation date
-regional_date <- as.Date("2022-01-06")
+regional_date <- as.Date("2022-02-12")
 
 # load case data
 daily <- load_local_data(regional_date) %>%
@@ -217,13 +217,14 @@ ggplot(piv_growth) +
 
 ```{r}
 posterior_summary %>%
+  as_tibble() %>%
   select(-rhat, -ess_bulk, -ess_tail) %>%
   filter(
     variable %in% c("gt_mean", "gt_sd", "ta", "voc_gt_mean", "voc_gt_sd")
   ) %>%
   mutate(across(where(is.numeric), round, digits = 1)) %>%
-  mutate(variant = case_when(variable %in% c("gt_mean", "gt_sd") ~ "Delta",
-                             TRUE ~ "Omicron")) %>%
+  mutate(variant = case_when(variable %in% c("gt_mean", "gt_sd") ~ "Omicron BA.1",
+                             TRUE ~ "Omicron BA.2")) %>%
   mutate(stratification = str_to_sentence(stratification),
          gt_type = str_to_sentence(gt_type)) %>%
   select(data_stratification = stratification,

writeup/summary.Rmd

@@ -60,11 +60,8 @@ pop <- load_population()
 sgtf <- load_results(date_latest)
 bias <- load_results(date_latest, type = "bias")
 
-sgtfposterior <- sgtf$posterior[,
-   variant_relationship := stringr::str_to_sentence(variant_relationship)
-]
 sgtf_posterior <- sgtf$posterior[
-  variant_relationship == "Correlated"
+  variant_relationship == "correlated"
 ]
 
 bias_posterior <- bias$posterior[variant_relationship == "correlated"]