Skip to content

Instantly share code, notes, and snippets.

@seabbs
Created March 9, 2021 11:18
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
Star You must be signed in to star a gist
Save seabbs/e36b6a8d3379783b4913f33a21cba7c9 to your computer and use it in GitHub Desktop.
Example script for processing samples.
# 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")))
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment