Example script for processing samples.

process-rt.R

# Packages ----------------------------------------------------------------
require(data.table, quietly = TRUE)
require(EpiNow2, quietly = TRUE)
require(purrr)
require(ggplot2)

# Target date -------------------------------------------------------------
creation_date <- Sys.Date()
extraction_date <- creation_date

# Get national -------------------------------------------------------------
deaths <- EpiNow2::get_regional_results(results_dir = here::here("forecast", "deaths", 
                                                                 "regional"), 
                                        date = extraction_date)$estimates$samples
linelist <- EpiNow2::get_regional_results(results_dir = here::here("forecast", "linelist-cases", 
                                                                   "regional"),
                                          date = extraction_date)$estimates$samples
admissions <- EpiNow2::get_regional_results(results_dir = here::here("forecast", "admissions", 
                                                                     "regional"), 
                                            date = extraction_date)$estimates$samples

combined <- data.table::rbindlist(list(deaths[, model := "EpiNow Deaths"],
                                       linelist[, model := "EpiNow Cases"],
                                       admissions[, model := "EpiNow Admissions"])) 

# Filter out forecasts ----------------------------------------------------
growth_measures <- data.table::copy(combined)[date <= creation_date][variable %in% c("R", "growth_rate", "infections")]

growth_measures[variable %in% "infections", variable := "incidence"]
infections <- data.table::copy(growth_measures)[variable %in% "incidence" & model == "EpiNow Admissions"][,
                                                type := "estimate"]
growth_measures <- data.table::rbindlist(list(growth_measures[!variable %in% "incidence"], infections))

# Work out doubling times -------------------------------------------------
doubling_time <- data.table::copy(growth_measures)[variable %in% "growth_rate", 
                                                    `:=`(variable = "doubling_time",
                                                         value = log(2) / value)]

# Convert gt sd -> variance -----------------------------------------------
kappa <- data.table::copy(combined)[variable == "gt_sd", 
                                    .(region, model, variable = "kappa", value = value ^ 2,
                                      date = list(unique(growth_measures[type %in% "estimate"]$date)))]

generation_time <- data.table::copy(combined)[variable == "gt_mean",
                                          .(region, model, variable = "mean_generation_time",
                                            value,  date = list(unique(growth_measures[type %in% "estimate"]$date)))]


gt_measures <- data.table::rbindlist(list(kappa, generation_time))
gt_measures <- gt_measures[, .(date = lubridate::as_date(unlist(date))), 
                              by = .(region, model, value, variable)][order(region, date)]
gt_measures <- gt_measures[!is.na(date)][, type := "estimate"]

# Extract measures of interest --------------------------------------------
combined <- data.table::rbindlist(
  list(growth_measures, doubling_time, gt_measures),
  fill = TRUE
)
combined[, c("parameter", "time", "sample", "strat") := NULL]

# Format as required ------------------------------------------------------
combined  <- combined[,
         .(Value = median(value, na.rm = TRUE), 
           `Quantile 0.05` = quantile(value, 0.05, na.rm = TRUE),
           `Quantile 0.1` = quantile(value, 0.1, na.rm = TRUE),
           `Quantile 0.15` = quantile(value, 0.15, na.rm = TRUE),
           `Quantile 0.2` = quantile(value, 0.2, na.rm = TRUE),
           `Quantile 0.25` = quantile(value, 0.25, na.rm = TRUE),
           `Quantile 0.3` = quantile(value, 0.3, na.rm = TRUE),
           `Quantile 0.35` = quantile(value, 0.35, na.rm = TRUE),
           `Quantile 0.4` = quantile(value, 0.4, na.rm = TRUE),
           `Quantile 0.45` = quantile(value, 0.45, na.rm = TRUE),
           `Quantile 0.5` = quantile(value, 0.5, na.rm = TRUE),
           `Quantile 0.55` = quantile(value, 0.55, na.rm = TRUE),
           `Quantile 0.6` = quantile(value, 0.6, na.rm = TRUE),
           `Quantile 0.65` = quantile(value, 0.65, na.rm = TRUE),
           `Quantile 0.7` = quantile(value, 0.7, na.rm = TRUE),
           `Quantile 0.75` = quantile(value, 0.75, na.rm = TRUE),
           `Quantile 0.8` = quantile(value, 0.8, na.rm = TRUE),
           `Quantile 0.85` = quantile(value, 0.85, na.rm = TRUE),
           `Quantile 0.9` = quantile(value, 0.9, na.rm = TRUE),
           `Quantile 0.95` = quantile(value, 0.95, na.rm = TRUE)), 
         by = c("region", "date", "variable", "model", "type")][,
            `:=`(
              Group = "LSHTM",
              `Creation Day` = lubridate::day(creation_date),
              `Creation Month` = lubridate::month(creation_date),
              `Creation Year` = lubridate::year(creation_date),
              `Day of Value` = lubridate::day(date),
              `Month of Value` = lubridate::month(date),
              `Year of Value` = lubridate::year(date),
              Geography = region,
              ValueType = variable,
              Model = model,
              Scenario = "Nowcast",
              ModelType = ifelse(model %in% "EpiNow Cases", "Cases",
                                 ifelse(model %in% "EpiNow Deaths", 
                                        "Deaths", "Cases")),
              Version = "2.0")
         ][, region := NULL][, variable := NULL][, model := NULL]

data.table::setcolorder(combined, c("Group", "Creation Day", "Creation Month", "Creation Year",
                              "Day of Value", "Month of Value", "Year of Value",
                              "Geography", "ValueType", "Model", "Scenario", "ModelType",  "Version"))
 
# Plot --------------------------------------------------------------------
last_estimate <-
  data.table::copy(combined)[type %in% "estimate"][,
                          .SD[date == max(date)], by = .(Geography,  Model, ValueType)]

## Linerange of just the last data point
linerange <- 
  ggplot(last_estimate, aes(x = Model, y = Value, col = Model)) +
  geom_linerange(aes(ymin = `Quantile 0.05`, ymax = `Quantile 0.95`),
                 alpha = 0.4, size = 5) +
  geom_linerange(aes(ymin = `Quantile 0.25`, ymax = `Quantile 0.75`),
                 alpha = 0.4, size = 5) +
  facet_wrap(ValueType ~ Geography, scales = "free_y") + 
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2") +
  cowplot::theme_cowplot() +
  theme(legend.position = "none") +
  labs(y = "Value", x = "Date", col = "Data")

ggsave(here::here("format-rt", "figures", "latest.png"),
       linerange, dpi = 330, height = 24, width = 24)


plot_timeseries <- function(combined, var = "R") {
  
  combined <- combined[, Model := ifelse(
                           Model %in% "EpiNow Cases", 
                           "Test-positive cases",
                           ifelse(Model %in% "EpiNow Admissions",
                                  "Hospital admissions", 
                                  "Deaths"))]
  ## Complete trend
  timeseries <- 
    ggplot(combined[type %in% "estimate"][ValueType %in% var], 
           aes(x = date, y = Value, col = Model, fill = Model)) +
    geom_ribbon(aes(ymin = `Quantile 0.05`, ymax = `Quantile 0.95`),
                alpha = 0.1, size = 0.1) +
    geom_ribbon(aes(col = NULL, ymin = `Quantile 0.25`, ymax = `Quantile 0.75`),
                alpha = 0.2) +
    geom_hline(yintercept = 1, linetype = 2) +
    facet_wrap(~ Geography, scales = "free_y", ncol = 2) + 
    scale_color_brewer(palette = "Dark2") +
    scale_fill_brewer(palette = "Dark2") +
    cowplot::theme_cowplot() +
    ggplot2::scale_x_date(date_breaks = "1 week", date_labels = "%b %d") +
    ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 90)) +
    theme(legend.position = "bottom") +
    labs(y = "R", x = "Date", col = "Data", fill = "Data")
  return(timeseries)
}

national_timeseries <- 
  plot_timeseries(combined[Geography %in% c("England", "Wales", "Scotland", "United Kingdom")])

ggsave(here::here("format-rt", "figures", "national-time-series.png"),
       national_timeseries, dpi = 330, height = 9, width = 9)

regional_timeseries <- 
  plot_timeseries(combined[!Geography %in% c("United Kingdom", "Northern Ireland")])

ggsave(here::here("format-rt", "figures", "regional-time-series.png"),
       regional_timeseries, dpi = 330, height = 12, width = 9)

# Report ------------------------------------------------------------------
## All estimates
data.table::fwrite(combined[Model %in% c("EpiNow Admissions",
                                         "EpiNow Cases")][, c("date", "type") := NULL],
                   here::here("format-rt", "data", "time-series", paste0(creation_date, "-time-series-r-lshtm.csv")))

## Estimates fully supported by data
data.table::fwrite(combined[Model %in% c("EpiNow Admissions",
                                         "EpiNow Cases")][type %in% "estimate"][, c("date", "type") := NULL],
                   here::here("format-rt", "data", "all", paste0(creation_date, "-r-lshtm.csv")))

# Report latest estimate only for each region -----------------------------

## Make Rt estimate for the UK be the Rt estimate for England
latest_r_estimate <- data.table::copy(combined)[ValueType %in% "R"][date == creation_date]
latest_r_estimate <- latest_r_estimate[Model == "EpiNow Cases"]
latest_r_estimate <- latest_r_estimate[, `Generation Process` := "Bayesian model based on the Cori et al. method but fit to latent 
                                       infections with a gaussian process to control temporal variation in the reproduction number."]
latest_r_estimate <- latest_r_estimate[, c("Group", "Generation Process",	"Creation Day",
                                           "Creation Month",	"Creation Year",	"Day of Value",
                                           "Month of Value",	"Year of Value",	"Geography",
                                           "ValueType",	"Quantile 0.05",	"Quantile 0.25", "Quantile 0.5",
                                           "Quantile 0.75",	"Quantile 0.95")]


data.table::fwrite(latest_r_estimate,
                   here::here("format-rt", "data", "current-r", paste0(creation_date, "-current-r-lshtm.csv")))