plots, results

compare/data-features.R

@@ -0,0 +1,53 @@
+# Raw case count features
+source("data/get-uk-data.R")
+
+data <- as.data.frame(data)
+
+data_regional <- split(data, data$region)
+
+# Mean difference between counts
+data_diffs <- purrr::map(data_regional, ~ 
+                          dplyr::mutate(., 
+                                        diff_case_death = .$cases_blend - .$deaths_blend,
+                                        diff_case_hosp = .$cases_blend - .$cases_hosp,
+                                        diff_hosp_death = .$cases_hosp - .$deaths_blend,
+                                        index = seq_len(nrow(.))) %>%
+                          dplyr::filter(date <= as.Date("2020-08-10")))
+
+data_mean <- data_diffs %>%
+  purrr::keep(names(.) %in% region_names$nhsregions) %>%
+  purrr::transpose() %>%
+  purrr::keep(names(.) %in% c("diff_case_death", "diff_case_hosp", "diff_hosp_death")) %>%
+  purrr::map_depth(.depth = 2, t.test, na.rm=T) %>%
+  purrr::map_depth(.depth = 2, ~ purrr::keep(., names(.) %in% c("estimate"))) %>%
+  purrr::map_depth(.depth = 2, ~ dplyr::bind_rows(.)) %>%
+  dplyr::bind_rows(., .id = "id")
+
+data_mean_min <- tidyr::pivot_longer(data_mean, cols = -id, values_to = "mean") %>%
+  dplyr::group_by(id) %>%
+  dplyr::filter(mean == min(mean))
+
+data_mean_max <- tidyr::pivot_longer(data_mean, cols = -id, values_to = "mean") %>%
+  dplyr::group_by(id) %>%
+  dplyr::filter(mean == max(mean))
+
+# Waves
+
+source("compare/wave-features.R")
+
+data_wave <- list(
+  cases = purrr::map(data_regional, ~ wave_features(.x, "cases_blend", 7)),
+  hosp = purrr::map(data_regional, ~ wave_features(.x, "cases_hosp", 7)),
+  deaths = purrr::map(data_regional, ~ wave_features(.x, "deaths_blend", 7)))
+
+data_wave <- data_wave %>%
+  purrr::transpose() %>%
+  purrr::discard(.p = names(.) %in% c("Scotland", "Wales", "Northern Ireland"))
+
+
+
+# Test positivity ---------------------------------------------------------
+
+t.test(pos_tests[pos_tests$region == "Midlands", "pos_perc"])
+
+

compare/plot-all.R

@@ -10,7 +10,7 @@ source("utils/utils.R")
 # Set global variables
 # consistent date axis:
 
-date_min <- as.Date("2020-04-13")
+date_min <- as.Date("2020-04-04")
 date_max <-as.Date("2020-08-20")
 theme_set(theme_classic(base_size = 12))
 

compare/plot-data.R

@@ -1,6 +1,8 @@
 require(dplyr)
 require(ggplot2)
 
+source("data/age-settings-data.R")
+
 # Get count data ----------------------------------------------------------
 # If data won't download, read in a saved hard copy of cleaned data - 23 August 2020
 if(class(try(source(here::here("data", "get-uk-data.R")))) == "try-error") {
@@ -45,15 +47,30 @@ data_ma <- standardised_data %>%
   dplyr::mutate('Data source' = dplyr::recode_factor(variable, 
                                                      "cases_blend" = "Cases",
                                                      "cases_hosp" = "Hospital admissions",
-                                                     "deaths_blend" = "Deaths"))
-
+                                                     "deaths_blend" = "Deaths")) %>%
+  # Add positivity rates (weekly)
+  dplyr::left_join(pos_tests, by = c("date", "region")) %>%
+  dplyr::mutate(pos_perc = factor(ifelse(!`Data source` == "Cases", NA, 
+                                         ifelse(pos_perc < 4.945 |
+                                                is.na(pos_perc), NA, ma))),
+                week = lubridate::week(date)) %>%
+  dplyr::group_by(region, week, `Data source`) %>%
+  tidyr::fill(pos_perc, .direction = "updown") %>%
+  dplyr::mutate(pos_perc = ifelse(is.na(pos_perc), NA, ma),
+                region = factor(region, levels = region_names$region_factor))
 
 # Plot --------------------------------------------------------------------
 
 # Plot 7 day MA for use with 2 rows
 plot_ma_only <- data_ma %>%
   ggplot() +
-  geom_line(aes(x = date, y = as.numeric(ma), colour = `Data source`)) +
+  geom_line(aes(x = date, y = as.numeric(ma), 
+                colour = `Data source`)) +
+  geom_point(aes(x = date, y = as.numeric(pos_perc),
+                 colour = `Data source`), 
+             shape = 3, size=0.9) +
+  scale_shape_discrete(solid=FALSE) +
+  geom_vline(xintercept = as.Date("2020-05-03"), lty = 3, colour = colours["Cases"]) +
   facet_wrap("region", nrow = 1, scales = "free_y") +
   cowplot::theme_cowplot() +
   coord_cartesian(xlim = c(date_min, date_max)) +
@@ -63,14 +80,20 @@ plot_ma_only <- data_ma %>%
         panel.spacing.y = unit(0.1, "cm"),
         axis.text.x = element_blank()) +
   guides(colour = FALSE) +
-  labs(y = "7-day MA", x = "")
+  labs(y = "7-day MA", x = NULL)
 
 
 # National - Plot 7 day MA for use with 2 rows
 plot_ma_only_national <- data_ma %>%
   dplyr::filter(region %in% c("England")) %>%
   ggplot() +
   geom_line(aes(x = date, y = as.numeric(ma), colour = `Data source`)) +
+  geom_point(aes(x = date, y = as.numeric(pos_perc),
+                 colour = `Data source`), 
+             shape = 3, size=0.9) +
+  geom_vline(xintercept = as.Date("2020-05-03"), 
+             lty = 3, colour = colours["Cases"],
+             alpha = 1) +
   cowplot::theme_cowplot() +
   coord_cartesian(xlim = c(date_min, date_max)) +
   scale_color_manual(values = colours) +
@@ -80,7 +103,7 @@ plot_ma_only_national <- data_ma %>%
         #axis.text.x = element_blank()
         ) +
   guides(colour = FALSE) +
-  labs(y = "7-day MA", x = "")
+  labs(y = "7-day MA", x = NULL)
 
 
 

compare/plot-ratios.R

@@ -49,7 +49,7 @@ plot_ratio_caseb_deathb <- summary_ratios %>%
   theme(strip.text.x = element_text(size = 18)) + ## Removes facet region name
   theme(axis.text.x = element_blank(),
         axis.text.y = element_text(size = 20)) + 
-  labs(y = "Rt(community) / Rt(deaths)", x = "")
+  labs(y = "Rt(community) / Rt(deaths)", x = NULL)
 
 # Hospital admissions on deaths
 plot_ratio_hosp_deathb <- summary_ratios %>%
@@ -72,7 +72,7 @@ plot_ratio_hosp_deathb <- summary_ratios %>%
   theme(strip.text.x = element_blank()) + 
   theme(axis.text.x = element_blank(),
         axis.text.y = element_text(size = 20)) + 
-  labs(y = "Rt(hospital) / Rt(deaths)",  x = "")
+  labs(y = "Rt(hospital) / Rt(deaths)",  x = NULL)
 
 # Cases by report date on hospital admissions
 plot_ratio_caseb_hosp <- summary_ratios %>%
@@ -95,7 +95,7 @@ plot_ratio_caseb_hosp <- summary_ratios %>%
         strip.text.x = element_blank(),
         axis.text.y = element_text(size = 20),
         axis.text.x = element_text(size = 20, angle = 45, hjust = 1)) + # Keep dates - bottom-most plot in grid
-  labs(y = "Rt(community) / Rt(hospital)", x = "")
+  labs(y = "Rt(community) / Rt(hospital)", x = NULL)
 
 
 
@@ -124,7 +124,7 @@ plot_national_ratio_caseb_deathb <- summary_ratios %>%
   scale_x_date(date_breaks = "1 month", date_labels = "%b") +
   cowplot::theme_cowplot() +
   theme(panel.spacing.x = unit(0.5, "cm")) +
-  labs(y = "", x = "")
+  labs(y = NULL, x = NULL)
 
 # Hospital admissions on deaths
 plot_national_hosp_deathb <- summary_ratios %>%
@@ -146,7 +146,7 @@ plot_national_hosp_deathb <- summary_ratios %>%
   scale_x_date(date_breaks = "1 month", date_labels = "%b") +
   cowplot::theme_cowplot() +
   theme(panel.spacing.x = unit(0.5, "cm")) +
-  labs(y = "",  x = "")
+  labs(y = NULL,  x = NULL)
 
 # Cases by report date on hospital admissions
 plot_national_caseb_hosp <- summary_ratios %>%
@@ -168,7 +168,7 @@ plot_national_caseb_hosp <- summary_ratios %>%
   scale_x_date(date_breaks = "1 month", date_labels = "%b") +
   cowplot::theme_cowplot() +
   theme(panel.spacing.x = unit(0.5, "cm")) + # Keep dates - bottom-most plot in grid
-  labs(y = "", x = "")
+  labs(y = NULL, x = NULL)
 
 
 

compare/plot-rt.R

@@ -53,7 +53,7 @@ plot_rt_only <- summary %>%
   theme(strip.text.x = element_blank()) +
   #theme(axis.text.x = element_blank()) +
   scale_x_date(date_breaks = "1 month", date_labels = "%b") +
-  labs(y = "R", x = "", col = "Data source", fill = "Data source") +
+  labs(y = "R", x = NULL, col = "Data source", fill = "Data source") +
   theme(legend.position = "none") +
   guides(fill = guide_legend(override.aes = list(alpha = 1)))
 
@@ -78,7 +78,7 @@ plot_rt_national <- summary %>%
   scale_fill_manual(values = colours) +
   scale_x_date(date_breaks = "1 month", date_labels = "%b") +
   cowplot::theme_cowplot() +
-  labs(y = "R", x = "", col = "Data source", fill = "Data source") +
+  labs(y = "R", x = NULL, col = "Data source", fill = "Data source") +
   theme(legend.position = "bottom") +
   guides(colour = guide_legend(override.aes = list(alpha = 1)))
 

compare/ratio-features.R

@@ -0,0 +1,31 @@
+# Ratio summary statistics
+summary_wide <- readRDS("rt-estimate/summary_wide.rds")
+
+library(dplyr)
+
+summary_ratio <- summary_wide
+#summary_ratio <- dplyr::filter(summary_ratio, date >= as.Date("2020-06-26"))
+summary_ratio <- split(summary_wide, summary_wide$region) 
+
+caseb_deathb_med <- summary_ratio %>%
+  purrr::map( ~ t.test(.x$caseb_deathb_med)) %>%
+  purrr::transpose() %>%
+  purrr::keep(., .p = names(.) %in% c("estimate", "conf.int")) %>%
+  purrr::transpose()
+
+caseb_hosp_med <- summary_ratio %>%
+  purrr::map( ~ t.test(.x$caseb_hosp_med)) %>%
+  purrr::transpose() %>%
+  purrr::keep(., .p = names(.) %in% c("estimate", "conf.int")) %>%
+  purrr::transpose()
+
+hosp_deathb_med <- summary_ratio %>%
+  purrr::map( ~ t.test(.x$hosp_deathb_med)) %>%
+  purrr::transpose() %>%
+  purrr::keep(., .p = names(.) %in% c("estimate", "conf.int")) %>%
+  purrr::transpose()
+
+
+median <- dplyr::select(summary_wide, date, region, dplyr::starts_with("median"))
+lower90 <- dplyr::select(summary_wide, date, region, dplyr::starts_with("lower_90"))
+upper90 <- dplyr::select(summary_wide, date, region, dplyr::starts_with("upper_90"))