pathways_paper_explore.Rmd

---
title: "Patient pathways paper - results"
author: "Quentin J. Leclerc, Gwenan M. Knight"
date: "`r format(Sys.time(), '%A %d %B %Y')`"
output:
  bookdown::pdf_document2:
    fig_width: 7.5
    fig_height: 4
---


# Setup 

## Load packages

```{r setup, include = F}
knitr::opts_chunk$set(echo = FALSE, message = F, warning = F,
                      fig.align = "center",
                      dpi = 900,
                      fig.path = "outputs/figures/",
                      dev = "pdf")

library(openxlsx)
library(dplyr)
library(linelist)
library(knitr)
library(kableExtra)
library(ggplot2)
library(reshape2)
library(ggpubr)
library(here)
library(cowplot)
library(patchwork)

source(here::here("code", "core_model_functions.R"))
source(here::here("code", "paths_table_to_blocks.R"))

set.seed(3210)

```

## Load England hospital data

```{r load_hospital_data}

#latest version of hosp data can be generated by sourcing the corresponding code:
#(generates a ton of objects, recommend clearing global envir after!)
#source(here::here("code", "import_hosp_data.R"))

hospital_data = read.csv(here::here("data", "uk_hospital_data.csv")) %>%
  mutate(date = as.Date(date))

cov_curve = hospital_data %>%
  filter(geography == "England") %>%
  filter(value_type == "hospital_inc")

eng_hospital_data = hospital_data %>%
  filter(geography == "England")

```

## Load UCH data

```{r load_uclh_data}

uclh_data = read.csv(here::here("data", "uclh_occupancy.csv")) %>%
  mutate(date = as.Date(date))

#need to extend cov_curve to fill in gaps and match max data date!
#in a future update, better to implement date column recognition in the model...
cov_curve_uclh = uclh_data %>%
  filter(value_type == "admission") %>%
  complete(date = seq.Date(min(date), max(uclh_data$date), by="day"), value_type) %>%
  mutate(value = replace_na(value, 0))


```

## Load model results

```{r load_model_results}

all_res = list.files(here::here("outputs", "results"))

for(filename in all_res){
  
  assign(strsplit(filename, split = ".rds")[[1]],
         readRDS(here::here("outputs", "results", filename)))
  
}

```

## Load NHS Trust results

```{r load_nhs_trusts}

best_fit_LoS_nhstrusts = read.csv(here::here("outputs", "tables", "best_fit_LoS_nhstrusts.csv"))

trust_to_region = read.csv(here::here("data", "trust_to_region.csv")) %>%
  select(trustname, nhs_region)

best_fit_LoS_nhstrusts$trustname = clean_variables(best_fit_LoS_nhstrusts)$region
best_fit_LoS_nhstrusts = merge(best_fit_LoS_nhstrusts, trust_to_region, by="trustname")

best_fit_LoS_nhstrusts = best_fit_LoS_nhstrusts %>%
  melt(id.vars = c("region", "trustname", "nhs_region"))

best_fit_LoS_nhstrusts = best_fit_LoS_nhstrusts %>%
  mutate(nhs_region = gsub("_", " ", nhs_region)) %>%
  mutate(nhs_region = stringr::str_to_title(nhs_region)) %>%
  mutate(nhs_region = gsub("And", "and", nhs_region)) %>%
  mutate(nhs_region = gsub("Of", "of", nhs_region))

best_fit_LoS_nhstrusts$variable = factor(best_fit_LoS_nhstrusts$variable, levels = c("CC", "Ward"))


```

# Visualise pathways

```{r pathsfig, fig.width = 15, fig.height = 6, dpi = 300}

data_orig = read.xlsx(here::here("data", "data_figs.xlsx"))

paths <- data_orig$path_name[1:7]
names(paths) <- data_orig$path[1:7]

stages <- c("Stage 1", "Stage 2", "Stage 3")
names(stages) <- c("1", "2", "3")

data <- subset(data_orig, path < 6)

data$path_f = factor(data$path, levels=c(2,1,4,5,3)) # To reorder to match proportions

p1 <- ggplot(data = subset(data, !is.na(bed) & source == "uclh"), aes(x=1, y = mean, fill = factor(bed)),na.rm = TRUE) + 
  coord_flip() + 
  geom_bar(stat="identity") + 
  geom_errorbar(aes(ymin = mean - sd, ymax = mean + sd)) + 
  facet_grid(path_f~step, labeller = labeller(path_f = paths, step = stages),  switch = "y", scales = "free") +  
  scale_x_continuous(labels = NULL, "") + 
  scale_y_continuous("Length of stay", limits = c(0,22)) + 
  theme_bw(base_size = 11)  + 
  scale_fill_discrete("Bed type", labels = c("CC","Ward"))+ 
  theme(axis.ticks.y = element_blank()) + 
  ggtitle("UCH")+
  theme(strip.text.y.left = element_text(angle = 0))

p2 <- ggplot(data = subset(data, !is.na(bed) & source == "cocin"), aes(x=1, y = mean, fill = factor(bed)),na.rm = TRUE) + 
  coord_flip() + 
  geom_bar(stat="identity") + 
  geom_errorbar(aes(ymin = pmax(0,mean - sd), ymax = pmin(mean + sd,22))) + 
  facet_grid(path_f~step, labeller = labeller(path_f = paths, step = stages),  switch = "y") + 
  scale_x_continuous(labels = NULL, "") + 
  scale_y_continuous("Length of stay", limits = c(0,22)) + 
  theme_bw(base_size = 11)  + 
  scale_fill_discrete("Bed type", labels = c("CC","Ward"))+ 
  theme(axis.ticks.y = element_blank()) + 
  ggtitle("COCIN") +
  theme(strip.text.y.left = element_text(angle = 0))

### Stacked bar chart of proportions
cbPalette <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")

data_prop <- data %>% filter(!is.na(proportion)) 

### Ordered by size
data_prop$path_name = with(data_prop, reorder(path_name, proportion))
data_prop$path_name = factor(data_prop$path_name)

p3 <- ggplot(data_prop,aes(x=source,y=proportion, fill = path_name)) + 
  geom_bar(position = "stack", stat = "identity") + 
  scale_fill_manual("Bed pathway",values=cbPalette) + 
  scale_x_discrete("Dataset", labels = c("COCIN","UCH")) + 
  scale_y_continuous("Proportion of admissions") +
  theme_bw(base_size = 11)

p3 + p1/p2 + plot_layout(ncol = 2, widths = c(1, 3)) + plot_annotation(tag_levels = 'A')




```

# Visualise model results

```{r uclh_res_melted}

Mres_average_LoS_uclh = res_average_LoS_uclh %>%
  select(-time) %>%
  melt(id.vars = "date")
Mres_average_LoS_uclh$grouping = "Model - average LoS"

Mres_full_pathways_uclh = res_full_pathways_uclh %>%
  select(-time) %>%
  melt(id.vars = "date")
Mres_full_pathways_uclh$grouping = "Model - full pathways"

Mres_all_uclh = rbind(Mres_average_LoS_uclh, Mres_full_pathways_uclh)

#separate mean and sd
results_mean_uclh = Mres_all_uclh %>%
  filter(!grepl("sd", variable)) %>%
  arrange(variable) %>%
  mutate(date = as.Date(date))

results_sd_uclh = Mres_all_uclh %>% 
  filter(grepl("sd", variable)) %>%
  arrange(variable) %>%
  mutate(date = as.Date(date))

results_sd_uclh$min = results_mean_uclh$value - 1.96*results_sd_uclh$value
results_sd_uclh$max = results_mean_uclh$value + 1.96*results_sd_uclh$value


```

```{r pointsinCIUCH}
#This tells us the percentage of data points that lie within the 95% CI of model

ward_points = results_sd_uclh %>%
  filter(variable == "Ward_sd",
         grouping == "Model - full pathways") %>%
  select(-value) %>%
  right_join({
    uclh_data %>%
      filter(value_type == "Ward")
  }, by = "date") %>%
  mutate(within = value <= max & value >= min)

round(table(ward_points$within)/length(ward_points$within), 2)

CC_points = results_sd_uclh %>%
  filter(variable == "CC_sd",
         grouping == "Model - full pathways") %>%
  select(-value) %>%
  right_join({
    uclh_data %>%
      filter(value_type == "CC")
  }, by = "date") %>%
  mutate(within = value <= max & value >= min)

round(table(CC_points$within)/length(CC_points$within), 2)


```

# Best-fit LoS distribution by NHS Trust

```{r losnhstrusts, fig.width = 9, fig.height = 6, dpi = 300}


#london midlands 
p1 = ggplot(data = subset(best_fit_LoS_nhstrusts, nhs_region %in% c("London", "Midlands")),
            aes(y = variable, x = value, colour = variable)) +
  geom_boxplot() +
  coord_flip() +
  theme_bw() +
  labs(x = "Length of stay (days)", y = "Bed type") +
  facet_wrap(~nhs_region) +
  guides(colour = F)


#north east north west
p2 = ggplot(data = subset(best_fit_LoS_nhstrusts, nhs_region %in% c("North East and Yorkshire", "North West")),
            aes(y = variable, x = value, colour = variable)) +
  geom_boxplot() +
  coord_flip() +
  theme_bw() +
  labs(x = "Length of stay (days)", y = "Bed type") +
  facet_wrap(~nhs_region) +
  guides(colour = F)

#south east east south west
p3 = ggplot(data = subset(best_fit_LoS_nhstrusts, nhs_region %in% c("South East", "South West", "East of England")),
            aes(y = variable, x = value, colour = variable)) +
  geom_boxplot() +
  coord_flip() +
  theme_bw() +
  labs(x = "Length of stay (days)", y = "Bed type") +
  facet_wrap(~nhs_region) +
  guides(colour = F)

#all together

plot_grid(p1,
          p2,
          p3,
          ncol = 1,
          rel_heights = c(1,1,1))

  

```

# Supplementary

```{r sqdiff_weibull}

sqdiff_weibull_average = results_mean_uclh %>%
  filter(grouping == "Model - average LoS") %>%
  select(-"grouping") %>%
  rename(value_type = variable) %>%
  mutate(dataset = "model") %>%
  rbind(cbind(uclh_data, dataset = "data")) %>%
  filter(value_type != "admission") %>%
  group_by(date, value_type) %>%
  summarise(sqdiff = diff(value)^2) %>%
  ungroup

sqdiff_weibull_full = results_mean_uclh %>%
  filter(grouping == "Model - full pathways") %>%
  select(-"grouping") %>%
  rename(value_type = variable) %>%
  mutate(dataset = "model") %>%
  rbind(cbind(uclh_data, dataset = "data")) %>%
  filter(value_type != "admission") %>%
  group_by(date, value_type) %>%
  summarise(sqdiff = diff(value)^2) %>%
  ungroup


```


```{r uclh_lognorm}

Mres_average_LoS_uclh = res_average_LoS_uclh_lognorm %>%
  select(-time) %>%
  melt(id.vars = "date")
Mres_average_LoS_uclh$grouping = "Model - average LoS"

Mres_full_pathways_uclh = res_full_pathways_uclh_lognorm %>%
  select(-time) %>%
  melt(id.vars = "date")
Mres_full_pathways_uclh$grouping = "Model - full pathways"

Mres_all_uclh = rbind(Mres_average_LoS_uclh, Mres_full_pathways_uclh)

#separate mean and sd
results_mean_uclh_lognorm = Mres_all_uclh %>%
  filter(!grepl("sd", variable)) %>%
  arrange(variable) %>%
  mutate(date = as.Date(date))

results_sd_uclh_lognorm = Mres_all_uclh %>% 
  filter(grepl("sd", variable)) %>%
  arrange(variable) %>%
  mutate(date = as.Date(date))

results_sd_uclh_lognorm$min = results_mean_uclh_lognorm$value - 1.96*results_sd_uclh_lognorm$value
results_sd_uclh_lognorm$max = results_mean_uclh_lognorm$value + 1.96*results_sd_uclh_lognorm$value

sqdiff_lognorm_average = results_mean_uclh_lognorm %>%
  filter(grouping == "Model - average LoS") %>%
  select(-"grouping") %>%
  rename(value_type = variable) %>%
  mutate(dataset = "model") %>%
  rbind(cbind(uclh_data, dataset = "data")) %>%
  filter(value_type != "admission") %>%
  group_by(date, value_type) %>%
  summarise(sqdiff = diff(value)^2) %>%
  ungroup

sqdiff_lognorm_full = results_mean_uclh_lognorm %>%
  filter(grouping == "Model - full pathways") %>%
  select(-"grouping") %>%
  rename(value_type = variable) %>%
  mutate(dataset = "model") %>%
  rbind(cbind(uclh_data, dataset = "data")) %>%
  filter(value_type != "admission") %>%
  group_by(date, value_type) %>%
  summarise(sqdiff = diff(value)^2) %>%
  ungroup


```

```{r uclh_weibull_round}

Mres_average_LoS_uclh = res_average_LoS_uclh_round %>%
  select(-time) %>%
  melt(id.vars = "date")
Mres_average_LoS_uclh$grouping = "Model - average LoS"

Mres_full_pathways_uclh = res_full_pathways_uclh_round %>%
  select(-time) %>%
  melt(id.vars = "date")
Mres_full_pathways_uclh$grouping = "Model - full pathways"

Mres_all_uclh = rbind(Mres_average_LoS_uclh, Mres_full_pathways_uclh)

#separate mean and sd
results_mean_uclh_round = Mres_all_uclh %>%
  filter(!grepl("sd", variable)) %>%
  arrange(variable) %>%
  mutate(date = as.Date(date))

results_sd_uclh_round = Mres_all_uclh %>% 
  filter(grepl("sd", variable)) %>%
  arrange(variable) %>%
  mutate(date = as.Date(date))

results_sd_uclh_round$min = results_mean_uclh_round$value - 1.96*results_sd_uclh_round$value
results_sd_uclh_round$max = results_mean_uclh_round$value + 1.96*results_sd_uclh_round$value

sqdiff_round_average = results_mean_uclh_round %>%
  filter(grouping == "Model - average LoS") %>%
  select(-"grouping") %>%
  rename(value_type = variable) %>%
  mutate(dataset = "model") %>%
  rbind(cbind(uclh_data, dataset = "data")) %>%
  filter(value_type != "admission") %>%
  group_by(date, value_type) %>%
  summarise(sqdiff = diff(value)^2) %>%
  ungroup

sqdiff_round_full = results_mean_uclh_round %>%
  filter(grouping == "Model - full pathways") %>%
  select(-"grouping") %>%
  rename(value_type = variable) %>%
  mutate(dataset = "model") %>%
  rbind(cbind(uclh_data, dataset = "data")) %>%
  filter(value_type != "admission") %>%
  group_by(date, value_type) %>%
  summarise(sqdiff = diff(value)^2) %>%
  ungroup


```

```{r uclh_lognorm_round}

Mres_average_LoS_uclh = res_average_LoS_uclh_lognorm_round %>%
  select(-time) %>%
  melt(id.vars = "date")
Mres_average_LoS_uclh$grouping = "Model - average LoS"

Mres_full_pathways_uclh = res_full_pathways_uclh_lognorm_round %>%
  select(-time) %>%
  melt(id.vars = "date")
Mres_full_pathways_uclh$grouping = "Model - full pathways"

Mres_all_uclh = rbind(Mres_average_LoS_uclh, Mres_full_pathways_uclh)

#separate mean and sd
results_mean_uclh_lognorm_round = Mres_all_uclh %>%
  filter(!grepl("sd", variable)) %>%
  arrange(variable) %>%
  mutate(date = as.Date(date))

results_sd_uclh_lognorm_round = Mres_all_uclh %>% 
  filter(grepl("sd", variable)) %>%
  arrange(variable) %>%
  mutate(date = as.Date(date))

results_sd_uclh_lognorm_round$min = results_mean_uclh_lognorm_round$value - 1.96*results_sd_uclh_lognorm_round$value
results_sd_uclh_lognorm_round$max = results_mean_uclh_lognorm_round$value + 1.96*results_sd_uclh_lognorm_round$value

sqdiff_lognorm_round_average = results_mean_uclh_lognorm_round %>%
  filter(grouping == "Model - average LoS") %>%
  select(-"grouping") %>%
  rename(value_type = variable) %>%
  mutate(dataset = "model") %>%
  rbind(cbind(uclh_data, dataset = "data")) %>%
  filter(value_type != "admission") %>%
  group_by(date, value_type) %>%
  summarise(sqdiff = diff(value)^2) %>%
  ungroup

sqdiff_lognorm_round_full = results_mean_uclh_lognorm_round %>%
  filter(grouping == "Model - full pathways") %>%
  select(-"grouping") %>%
  rename(value_type = variable) %>%
  mutate(dataset = "model") %>%
  rbind(cbind(uclh_data, dataset = "data")) %>%
  filter(value_type != "admission") %>%
  group_by(date, value_type) %>%
  summarise(sqdiff = diff(value)^2) %>%
  ungroup


```


```{r sqdifftable}

sqdiff_distributions = data.frame(rounding = c(rep("Round", 4), rep("Ceiling", 4)),
                                  distribution = rep(c("Weibull", "Weibull", "Lognormal", "Lognormal"),2),
                                  pathway = c("Average LoS", "Bed pathways",
                                              "Average LoS", "Bed pathways",
                                              "Average LoS", "Bed pathways",
                                              "Average LoS", "Bed pathways"),
                                  sqdiff = c(sum(sqdiff_round_average$sqdiff),
                                             sum(sqdiff_round_full$sqdiff),
                                             sum(sqdiff_lognorm_round_average$sqdiff),
                                             sum(sqdiff_lognorm_round_full$sqdiff),
                                             sum(sqdiff_weibull_average$sqdiff),
                                             sum(sqdiff_weibull_full$sqdiff),
                                             sum(sqdiff_lognorm_average$sqdiff),
                                             sum(sqdiff_lognorm_full$sqdiff)))

kable(sqdiff_distributions, format = "latex", booktabs = T,
      caption = "Squared difference between model results for average LoS or bed pathways, and Weibull or Lognormal 
      distribution.",
      col.names = c("Rounding",
                    "Distribution",
                    "LoS used",
                    "Squared difference")) %>%
  kable_styling(latex_options = c("striped"))

```

```{r lognormweibull, fig.width = 14, fig.height = 7, dpi = 300}

weibullpar = mixdist::weibullpar(3.98, 3.60, 0)
weibull_LoS = ceiling(rweibull(1000, shape = weibullpar$shape, scale = weibullpar$scale))

lognorm_LoS = ceiling(rlnorm(1000,
                           log(3.98),
                           log(3.60)))

#create plot C
plotC = ggplot() +
  geom_density(aes(weibull_LoS, fill = "Weibull"), alpha = 0.2) +
  geom_density(aes(lognorm_LoS, fill = "Lognormal"), alpha = 0.2) +
  scale_x_continuous(breaks = seq(1,31,4), limits = c(0,31)) +
  labs(x = "LoS (days)", y = "Density", fill = "Distribution:")+
  theme_bw()

#create plot A
plotA = ggplot() +
  geom_line(data = results_mean_uclh %>% filter(grouping == "Model - average LoS"),
            aes(date, value, colour = variable, linetype = "Model"), size = 1) +
  geom_ribbon(data = subset(results_sd_uclh, variable == "CC_sd") %>%
                filter(grouping == "Model - average LoS"), 
              aes(date,
                  ymin = min,
                  ymax = max,
                  fill = "CC"),
              alpha = 0.3) +
  geom_ribbon(data = subset(results_sd_uclh, variable == "Ward_sd")%>%
                filter(grouping == "Model - average LoS"),
              aes(date,
                  ymin = min,
                  ymax = max,
                  fill = "Ward"),
              alpha = 0.3) +
  geom_line(data = uclh_data %>% filter(value_type == "Ward"),
            aes(date, value, linetype = "Data", colour = "Ward"), size = 1) +
  geom_line(data = uclh_data %>% filter(value_type == "CC"),
            aes(date, value, linetype = "Data", colour = "CC"), size = 1) +
  labs(colour = "Bed type:", linetype = "Value:",
       x = "Time (days)", y = "Bed occupancy") +
  guides(fill = FALSE) +
  scale_x_date(breaks = "1 week", date_labels = format("%d %b %Y")) +
  scale_y_continuous(limits = c(-1.5,70)) +
  theme_bw()

#create plot B
plotB = ggplot() +
  geom_line(data = results_mean_uclh_lognorm %>% filter(grouping == "Model - average LoS"),
            aes(date, value, colour = variable, linetype = "Model"), size = 1) +
  geom_ribbon(data = subset(results_sd_uclh_lognorm, variable == "CC_sd") %>%
                filter(grouping == "Model - average LoS"), 
              aes(date,
                  ymin = min,
                  ymax = max,
                  fill = "CC"),
              alpha = 0.3) +
  geom_ribbon(data = subset(results_sd_uclh_lognorm, variable == "Ward_sd")%>%
                filter(grouping == "Model - average LoS"),
              aes(date,
                  ymin = min,
                  ymax = max,
                  fill = "Ward"),
              alpha = 0.3) +
  geom_line(data = uclh_data %>% filter(value_type == "Ward"),
            aes(date, value, linetype = "Data", colour = "Ward"), size = 1) +
  geom_line(data = uclh_data %>% filter(value_type == "CC"),
            aes(date, value, linetype = "Data", colour = "CC"), size = 1) +
  labs(colour = "Bed type:", linetype = "Value:",
       x = "Time (days)", y = "Bed occupancy") +
  guides(fill = FALSE) +
  scale_x_date(breaks = "1 week", date_labels = format("%d %b %Y")) +
  scale_y_continuous(limits = c(-1.5,70)) +
  theme_bw()

legend = get_legend(plotA)
plot_grid(
  plot_grid(
    plotA+theme(legend.position = "none"), 
    plotB+theme(legend.position = "none"), 
    legend, 
    labels = c("A)", "B)"), label_size = 12,
    hjust = -0.3,
    rel_widths = c(1, 1, 0.3), nrow = 1),
  plot_grid(plotC, 
            labels = "C)", label_size = 12,
            hjust = -0.3, 
            nrow = 1),
  ncol = 1)

```

# Look at COCIN data

## Plot all the simulations together

```{r plot_all_together}

Mres_average_LoS = res_average_LoS_eng %>%
  select(-time) %>%
  melt(id.vars = "date")
Mres_average_LoS$grouping = "Model - average LoS"

Mres_full_pathways = res_full_pathways_eng %>%
  select(-time) %>%
  melt(id.vars = "date")
Mres_full_pathways$grouping = "Model - full pathways"

Mres_all = rbind(Mres_average_LoS, Mres_full_pathways)

#separate mean and sd
results_mean = Mres_all %>%
  filter(!grepl("sd", variable)) %>%
  arrange(variable) %>%
  mutate(date = as.Date(date))

results_sd = Mres_all %>% 
  filter(grepl("sd", variable)) %>%
  arrange(variable) %>%
  mutate(date = as.Date(date))

results_sd$min = results_mean$value - 1.96*results_sd$value
results_sd$max = results_mean$value + 1.96*results_sd$value


#same but for rounded results

Mres_average_LoS = res_average_LoS_eng_round %>%
  select(-time) %>%
  melt(id.vars = "date")
Mres_average_LoS$grouping = "Model - average LoS"

Mres_full_pathways = res_full_pathways_eng_round %>%
  select(-time) %>%
  melt(id.vars = "date")
Mres_full_pathways$grouping = "Model - full pathways"

Mres_all = rbind(Mres_average_LoS, Mres_full_pathways)

results_mean_round = Mres_all %>%
  filter(!grepl("sd", variable)) %>%
  arrange(variable) %>%
  mutate(date = as.Date(date))

results_sd_round = Mres_all %>% 
  filter(grepl("sd", variable)) %>%
  arrange(variable) %>%
  mutate(date = as.Date(date))

results_sd_round$min = results_mean_round$value - 1.96*results_sd_round$value
results_sd_round$max = results_mean_round$value + 1.96*results_sd_round$value


```

## Massive figure for both datasets

```{r massiveplot, fig.width = 12, fig.height = 7, dpi = 300}

#create plot A
plotA = ggplot() +
  geom_line(data = results_mean_uclh %>% filter(grouping == "Model - average LoS"),
            aes(date, value, colour = variable, linetype = "Model"), size = 1) +
  geom_ribbon(data = subset(results_sd_uclh, variable == "CC_sd") %>%
                filter(grouping == "Model - average LoS"), 
              aes(date,
                  ymin = min,
                  ymax = max,
                  fill = "CC"),
              alpha = 0.3) +
  geom_ribbon(data = subset(results_sd_uclh, variable == "Ward_sd")%>%
                filter(grouping == "Model - average LoS"),
              aes(date,
                  ymin = min,
                  ymax = max,
                  fill = "Ward"),
              alpha = 0.3) +
  geom_line(data = uclh_data %>% filter(value_type == "Ward"),
            aes(date, value, linetype = "Data", colour = "Ward"), size = 1) +
  geom_line(data = uclh_data %>% filter(value_type == "CC"),
            aes(date, value, linetype = "Data", colour = "CC"), size = 1) +
  labs(colour = "Bed type:", linetype = "Value:",
       x = "Time (days)", y = "Bed occupancy") +
  guides(fill = FALSE) +
  scale_x_date(breaks = "1 week", date_labels = format("%d %b %Y")) +
  scale_y_continuous(limits = c(-1.5,65)) +
  theme_bw()

#create plot B
plotB = ggplot() +
  geom_line(data = results_mean_uclh %>% filter(grouping == "Model - full pathways"),
            aes(date, value, colour = variable, linetype = "Model"), size = 1) +
  geom_ribbon(data = subset(results_sd_uclh, variable == "CC_sd") %>%
                filter(grouping == "Model - full pathways"), 
              aes(date,
                  ymin = min,
                  ymax = max,
                  fill = "CC"),
              alpha = 0.3) +
  geom_ribbon(data = subset(results_sd_uclh, variable == "Ward_sd")%>%
                filter(grouping == "Model - full pathways"),
              aes(date,
                  ymin = min,
                  ymax = max,
                  fill = "Ward"),
              alpha = 0.3) +
  geom_line(data = uclh_data %>% filter(value_type == "Ward"),
            aes(date, value, linetype = "Data", colour = "Ward"), size = 1) +
  geom_line(data = uclh_data %>% filter(value_type == "CC"),
            aes(date, value, linetype = "Data", colour = "CC"), size = 1) +
  labs(colour = "Bed type:", linetype = "Value:",
       x = "Time (days)", y = "Bed occupancy") +
  guides(fill = FALSE) +
  scale_x_date(breaks = "1 week", date_labels = format("%d %b %Y")) +
  scale_y_continuous(limits = c(-1.5,65)) +
  theme_bw()

#create plot C
plotC = ggplot() +
  geom_line(data = results_mean %>% filter(grouping == "Model - average LoS"),
            aes(date, value, colour = variable, linetype = "Model"), size = 1) +
  geom_ribbon(data = subset(results_sd, variable == "CC_sd") %>%
                filter(grouping == "Model - average LoS"), 
              aes(date,
                  ymin = min,
                  ymax = max,
                  fill = "CC"),
              alpha = 0.3) +
  geom_ribbon(data = subset(results_sd, variable == "Ward_sd")%>%
                filter(grouping == "Model - average LoS"),
              aes(date,
                  ymin = min,
                  ymax = max,
                  fill = "Ward"),
              alpha = 0.3) +
  geom_line(data = eng_hospital_data %>% filter(value_type == "hospital_prev"),
            aes(date, value, linetype = "Data", colour = "Ward"), size = 1) +
  geom_line(data = eng_hospital_data %>% filter(value_type == "icu_prev"),
            aes(date, value, linetype = "Data", colour = "CC"), size = 1) +
  labs(colour = "Bed type:", linetype = "Value:",
       x = "Time (days)", y = "Bed occupancy") +
  guides(fill = FALSE) +
  scale_x_date(breaks = "1 month", date_labels = format("%d %b %Y")) +
  scale_y_continuous(limits = c(-1.5,21000)) +
  theme_bw()

#create plot D
plotD = ggplot() +
  geom_line(data = results_mean %>% filter(grouping == "Model - full pathways"),
            aes(date, value, colour = variable, linetype = "Model"), size = 1) +
  geom_ribbon(data = subset(results_sd, variable == "CC_sd") %>%
                filter(grouping == "Model - full pathways"), 
              aes(date,
                  ymin = min,
                  ymax = max,
                  fill = "CC"),
              alpha = 0.3) +
  geom_ribbon(data = subset(results_sd, variable == "Ward_sd")%>%
                filter(grouping == "Model - full pathways"),
              aes(date,
                  ymin = min,
                  ymax = max,
                  fill = "Ward"),
              alpha = 0.3) +
  geom_line(data = eng_hospital_data %>% filter(value_type == "hospital_prev"),
            aes(date, value, linetype = "Data", colour = "Ward"), size = 1) +
  geom_line(data = eng_hospital_data %>% filter(value_type == "icu_prev"),
            aes(date, value, linetype = "Data", colour = "CC"), size = 1) +
  labs(colour = "Bed type:", linetype = "Value:",
       x = "Time (days)", y = "Bed occupancy") +
  guides(fill = FALSE) +
  scale_x_date(breaks = "1 month", date_labels = format("%d %b %Y")) +
  scale_y_continuous(limits = c(-1.5,21000)) +
  theme_bw()


legend = get_legend(plotA + theme(legend.position = "bottom"))
plot_grid(
  NULL,
  plot_grid(
    plotA+theme(legend.position = "none"), plotB+theme(legend.position = "none"),
    NULL, NULL,
    plotC+theme(legend.position = "none"), plotD+theme(legend.position = "none"),
    labels = c("A) UCH, average LoS",
               "B) UCH, bed pathways",
               "", "",
               "C) COCIN, average LoS",
               "D) COCIN, bed pathways"), 
    label_size = 12, hjust = -0.3, vjust = 0,
    nrow = 3,
    rel_heights = c(1, 0.05, 1)
  ),
  legend,
  nrow = 3,
  rel_heights = c(0.05, 1, 0.1))


```

## Supplementary - Massive figure for both datasets, round instead of ceiling

```{r massiveplot_rounded, fig.width = 12, fig.height = 7, dpi = 300}

#create plot A
plotA = ggplot() +
  geom_line(data = results_mean_uclh_round %>% filter(grouping == "Model - average LoS"),
            aes(date, value, colour = variable, linetype = "Model"), size = 1) +
  geom_ribbon(data = subset(results_sd_uclh_round, variable == "CC_sd") %>%
                filter(grouping == "Model - average LoS"), 
              aes(date,
                  ymin = min,
                  ymax = max,
                  fill = "CC"),
              alpha = 0.3) +
  geom_ribbon(data = subset(results_sd_uclh_round, variable == "Ward_sd")%>%
                filter(grouping == "Model - average LoS"),
              aes(date,
                  ymin = min,
                  ymax = max,
                  fill = "Ward"),
              alpha = 0.3) +
  geom_line(data = uclh_data %>% filter(value_type == "Ward"),
            aes(date, value, linetype = "Data", colour = "Ward"), size = 1) +
  geom_line(data = uclh_data %>% filter(value_type == "CC"),
            aes(date, value, linetype = "Data", colour = "CC"), size = 1) +
  labs(colour = "Bed type:", linetype = "Value:",
       x = "Time (days)", y = "Bed occupancy") +
  guides(fill = FALSE) +
  scale_x_date(breaks = "1 week", date_labels = format("%d %b %Y")) +
  theme_bw()

#create plot B
plotB = ggplot() +
  geom_line(data = results_mean_uclh_round %>% filter(grouping == "Model - full pathways"),
            aes(date, value, colour = variable, linetype = "Model"), size = 1) +
  geom_ribbon(data = subset(results_sd_uclh_round, variable == "CC_sd") %>%
                filter(grouping == "Model - full pathways"), 
              aes(date,
                  ymin = min,
                  ymax = max,
                  fill = "CC"),
              alpha = 0.3) +
  geom_ribbon(data = subset(results_sd_uclh_round, variable == "Ward_sd")%>%
                filter(grouping == "Model - full pathways"),
              aes(date,
                  ymin = min,
                  ymax = max,
                  fill = "Ward"),
              alpha = 0.3) +
  geom_line(data = uclh_data %>% filter(value_type == "Ward"),
            aes(date, value, linetype = "Data", colour = "Ward"), size = 1) +
  geom_line(data = uclh_data %>% filter(value_type == "CC"),
            aes(date, value, linetype = "Data", colour = "CC"), size = 1) +
  labs(colour = "Bed type:", linetype = "Value:",
       x = "Time (days)", y = "Bed occupancy") +
  guides(fill = FALSE) +
  scale_x_date(breaks = "1 week", date_labels = format("%d %b %Y")) +
  theme_bw()

#create plot C
plotC = ggplot() +
  geom_line(data = results_mean_round %>% filter(grouping == "Model - average LoS"),
            aes(date, value, colour = variable, linetype = "Model"), size = 1) +
  geom_ribbon(data = subset(results_sd_round, variable == "CC_sd") %>%
                filter(grouping == "Model - average LoS"), 
              aes(date,
                  ymin = min,
                  ymax = max,
                  fill = "CC"),
              alpha = 0.3) +
  geom_ribbon(data = subset(results_sd_round, variable == "Ward_sd")%>%
                filter(grouping == "Model - average LoS"),
              aes(date,
                  ymin = min,
                  ymax = max,
                  fill = "Ward"),
              alpha = 0.3) +
  geom_line(data = eng_hospital_data %>% filter(value_type == "hospital_prev"),
            aes(date, value, linetype = "Data", colour = "Ward"), size = 1) +
  geom_line(data = eng_hospital_data %>% filter(value_type == "icu_prev"),
            aes(date, value, linetype = "Data", colour = "CC"), size = 1) +
  labs(colour = "Bed type:", linetype = "Value:",
       x = "Time (days)", y = "Bed occupancy") +
  guides(fill = FALSE) +
  scale_x_date(breaks = "1 month", date_labels = format("%d %b %Y")) +
  theme_bw()

#create plot D
plotD = ggplot() +
  geom_line(data = results_mean_round %>% filter(grouping == "Model - full pathways"),
            aes(date, value, colour = variable, linetype = "Model"), size = 1) +
  geom_ribbon(data = subset(results_sd_round, variable == "CC_sd") %>%
                filter(grouping == "Model - full pathways"), 
              aes(date,
                  ymin = min,
                  ymax = max,
                  fill = "CC"),
              alpha = 0.3) +
  geom_ribbon(data = subset(results_sd_round, variable == "Ward_sd")%>%
                filter(grouping == "Model - full pathways"),
              aes(date,
                  ymin = min,
                  ymax = max,
                  fill = "Ward"),
              alpha = 0.3) +
  geom_line(data = eng_hospital_data %>% filter(value_type == "hospital_prev"),
            aes(date, value, linetype = "Data", colour = "Ward"), size = 1) +
  geom_line(data = eng_hospital_data %>% filter(value_type == "icu_prev"),
            aes(date, value, linetype = "Data", colour = "CC"), size = 1) +
  labs(colour = "Bed type:", linetype = "Value:",
       x = "Time (days)", y = "Bed occupancy") +
  guides(fill = FALSE) +
  scale_x_date(breaks = "1 month", date_labels = format("%d %b %Y")) +
  theme_bw()


legend = get_legend(plotA + theme(legend.position = "bottom"))
plot_grid(
  plot_grid(
    plotA+theme(legend.position = "none"), plotB+theme(legend.position = "none"),
    plotC+theme(legend.position = "none"), plotD+theme(legend.position = "none"),
    labels = c("A)", "B)", "C)", "D)"), 
    label_size = 12, hjust = -0.3
  ),
  legend,
  nrow = 2,
  rel_heights = c(1, 0.1))


```


## Summary table of peaks

```{r tablesummary}

#uclh
uclh_data_ward = uclh_data %>%
  filter(value_type == "Ward")
uclh_data_CC = uclh_data %>%
  filter(value_type == "CC")

t_max_average_LoS_ward_uclh = median(res_average_LoS_uclh_peak$Ward_peak_time)
t_max_full_pathways_ward_uclh = median(res_full_pathways_uclh_peak$Ward_peak_time)
t_max_data_ward_uclh = which.max(uclh_data_ward$value)

t_max_average_LoS_CC_uclh = median(res_average_LoS_uclh_peak$CC_peak_time)
t_max_full_pathways_CC_uclh = median(res_full_pathways_uclh_peak$CC_peak_time)
t_max_data_CC_uclh = which.max(uclh_data_CC$value)


#cocin
eng_hospital_data_ward = eng_hospital_data %>%
  filter(value_type == "hospital_prev")
eng_hospital_data_CC = eng_hospital_data %>%
  filter(value_type == "icu_prev")

t_max_average_LoS_ward = median(res_average_LoS_eng_peak$Ward_peak_time)
t_max_full_pathways_ward = median(res_full_pathways_eng_peak$Ward_peak_time)
t_max_data_ward = which.max(eng_hospital_data_ward$value)

t_max_average_LoS_CC = median(res_average_LoS_eng_peak$CC_peak_time)
t_max_full_pathways_CC = median(res_full_pathways_eng_peak$CC_peak_time)
t_max_data_CC = which.max(eng_hospital_data_CC$value)


summary_results = data.frame(Dataset = c(rep("UCH", 3), rep("COCIN", 3)),
                             Grouping = rep(c("Average LoS",
                                              "Bed pathways",
                                              "Data"), 2),
                             ward_t_max = c(res_average_LoS_uclh$date[t_max_average_LoS_ward_uclh],
                                            res_full_pathways_uclh$date[t_max_full_pathways_ward_uclh],
                                            uclh_data_ward$date[t_max_data_ward_uclh],
                                            res_average_LoS_eng$date[t_max_average_LoS_ward],
                                            res_full_pathways_eng$date[t_max_full_pathways_ward],
                                            eng_hospital_data_ward$date[t_max_data_ward]),
                             ward_t_quant = c(res_average_LoS_uclh$date[quantile(res_average_LoS_uclh_peak$Ward_peak_time, 0.25)],
                                              res_full_pathways_uclh$date[quantile(res_full_pathways_uclh_peak$Ward_peak_time, 0.25)],
                                              NA,
                                              res_average_LoS_eng$date[quantile(res_average_LoS_eng_peak$Ward_peak_time, 0.25)],
                                              res_full_pathways_eng$date[quantile(res_full_pathways_eng_peak$Ward_peak_time, 0.25)],
                                              NA),
                             ward_t_quant2 = c(res_average_LoS_uclh$date[quantile(res_average_LoS_uclh_peak$Ward_peak_time, 0.75)],
                                               res_full_pathways_uclh$date[quantile(res_full_pathways_uclh_peak$Ward_peak_time, 0.75)],
                                               NA,
                                               res_average_LoS_eng$date[quantile(res_average_LoS_eng_peak$Ward_peak_time, 0.75)],
                                               res_full_pathways_eng$date[quantile(res_full_pathways_eng_peak$Ward_peak_time, 0.75)],
                                               NA),
                             
                             ward_max = c(mean(res_average_LoS_uclh_peak$Ward_peak),
                                          mean(res_full_pathways_uclh_peak$Ward_peak),
                                          uclh_data_ward$value[t_max_data_ward_uclh],
                                          mean(res_average_LoS_eng_peak$Ward_peak),
                                          mean(res_full_pathways_eng_peak$Ward_peak),
                                          eng_hospital_data_ward$value[t_max_data_ward]),
                             ward_max_sd = c(sd(res_average_LoS_uclh_peak$Ward_peak),
                                             sd(res_full_pathways_uclh_peak$Ward_peak),
                                             NA,
                                             sd(res_average_LoS_eng_peak$Ward_peak),
                                             sd(res_full_pathways_eng_peak$Ward_peak),
                                             NA),
                             
                             cc_t_max = c(res_average_LoS_uclh$date[t_max_average_LoS_CC_uclh],
                                          res_full_pathways_uclh$date[t_max_full_pathways_CC_uclh],
                                          uclh_data_CC$date[t_max_data_CC_uclh],
                                          res_average_LoS_eng$date[t_max_average_LoS_CC],
                                          res_full_pathways_eng$date[t_max_full_pathways_CC],
                                          eng_hospital_data_CC$date[t_max_data_CC]),
                             cc_t_quant = c(res_average_LoS_uclh$date[quantile(res_average_LoS_uclh_peak$CC_peak_time, 0.25)],
                                            res_full_pathways_uclh$date[quantile(res_full_pathways_uclh_peak$CC_peak_time, 0.25)],
                                            NA,
                                            res_average_LoS_eng$date[quantile(res_average_LoS_eng_peak$CC_peak_time, 0.25)],
                                            res_full_pathways_eng$date[quantile(res_full_pathways_eng_peak$CC_peak_time, 0.25)],
                                            NA),
                             cc_t_quant2 = c(res_average_LoS_uclh$date[quantile(res_average_LoS_uclh_peak$CC_peak_time, 0.75)],
                                             res_full_pathways_uclh$date[quantile(res_full_pathways_uclh_peak$CC_peak_time, 0.75)],
                                             NA,
                                             res_average_LoS_eng$date[quantile(res_average_LoS_eng_peak$CC_peak_time, 0.75)],
                                             res_full_pathways_eng$date[quantile(res_full_pathways_eng_peak$CC_peak_time, 0.75)],
                                             NA),
                             
                             cc_max = c(mean(res_average_LoS_uclh_peak$CC_peak),
                                        mean(res_full_pathways_uclh_peak$CC_peak),
                                        uclh_data_CC$value[t_max_data_CC_uclh],
                                        mean(res_average_LoS_eng_peak$CC_peak),
                                        mean(res_full_pathways_eng_peak$CC_peak),
                                        eng_hospital_data_CC$value[t_max_data_CC]),
                             cc_max_sd = c(sd(res_average_LoS_uclh_peak$CC_peak),
                                           sd(res_full_pathways_uclh_peak$CC_peak),
                                           NA,
                                           sd(res_average_LoS_eng_peak$CC_peak),
                                           sd(res_full_pathways_eng_peak$CC_peak),
                                           NA))

kable(summary_results, format = "latex", booktabs = T,
      caption = "Summary bed occupancy for the 2 models and the data. CC: critical care (ICU or HDU).",
      col.names = c("Dataset", "LoS used",
                    "Peak time", "Peak IQR1", "Peak IQR2", "Mean beds", "SD",
                    "Peak time", "Peak IQR1", "Peak IQR2", "Mean beds", "SD")) %>%
  add_header_above(c(" " = 2, "Ward peak bed occupancy" = 5, "CC peak bed occupancy" = 5)) %>%
  kable_styling(latex_options = c("striped", "scale_down"))


```

## By region

```{r region_plots_average, fig.width = 9, fig.height = 6, dpi = 300}

#london midlands 
p1 = ggplot() +
  geom_line(data = subset(res_regions_eng, region %in% c("London", "Midlands")),
            aes(date, Ward, colour = "Ward", linetype = "2Model"), size = 1) +
  geom_ribbon(data = subset(res_regions_eng, region %in% c("London", "Midlands")),
              aes(date,
                  ymin = Ward - 1.96*Ward_sd,
                  ymax = Ward + 1.96*Ward_sd,
                  fill = "Ward"), alpha = 0.3) +
  geom_line(data = subset(data_regions_eng, region %in% c("London", "Midlands")),
            aes(date, Ward, colour = "Ward", linetype = "1Data"), size = 1) +
  geom_line(data = subset(res_regions_eng, region %in% c("London", "Midlands")),
            aes(date, CC, colour = "CC", linetype = "2Model"), size = 1) +
  geom_ribbon(data = subset(res_regions_eng, region %in% c("London", "Midlands")),
              aes(date,
                  ymin = CC - 1.96*CC_sd,
                  ymax = CC + 1.96*CC_sd,
                  fill = "CC"), alpha = 0.3) +
  geom_line(data = subset(data_regions_eng, region %in% c("London", "Midlands")),
            aes(date, CC, colour = "CC", linetype = "1Data"), size = 1) +
  scale_linetype_discrete(labels = c("Data", "Model")) +
  facet_wrap(~region) +
  guides(fill = F) +
  labs(x = "Date", y = "Bed occupancy", color = "Bed type:", linetype = "Value:") +
  theme_bw()

#north east north west
p2 = ggplot() +
  geom_line(data = subset(res_regions_eng, region %in% c("North East and Yorkshire", "North West")),
            aes(date, Ward, colour = "Ward", linetype = "2Model"), size = 1) +
  geom_ribbon(data = subset(res_regions_eng, region %in% c("North East and Yorkshire", "North West")),
              aes(date,
                  ymin = Ward - 1.96*Ward_sd,
                  ymax = Ward + 1.96*Ward_sd,
                  fill = "Ward"), alpha = 0.3) +
  geom_line(data = subset(data_regions_eng, region %in% c("North East and Yorkshire", "North West")),
            aes(date, Ward, colour = "Ward", linetype = "1Data"), size = 1) +
  geom_line(data = subset(res_regions_eng, region %in% c("North East and Yorkshire", "North West")),
            aes(date, CC, colour = "CC", linetype = "2Model"), size = 1) +
  geom_ribbon(data = subset(res_regions_eng, region %in% c("North East and Yorkshire", "North West")),
              aes(date,
                  ymin = CC - 1.96*CC_sd,
                  ymax = CC + 1.96*CC_sd,
                  fill = "CC"), alpha = 0.3) +
  geom_line(data = subset(data_regions_eng, region %in% c("North East and Yorkshire", "North West")),
            aes(date, CC, colour = "CC", linetype = "1Data"), size = 1) +
  scale_linetype_discrete(labels = c("Data", "Model")) +
  facet_wrap(~region) +
  guides(fill = F) +
  labs(x = "Date", y = "Bed occupancy", color = "Bed type:", linetype = "Value:") +
  theme_bw()

#south east east south west
p3 = ggplot() +
  geom_line(data = subset(res_regions_eng, region %in% c("South East", "South West", "East of England")),
            aes(date, Ward, colour = "Ward", linetype = "2Model"), size = 1) +
  geom_ribbon(data = subset(res_regions_eng, region %in% c("South East", "South West", "East of England")),
              aes(date,
                  ymin = Ward - 1.96*Ward_sd,
                  ymax = Ward + 1.96*Ward_sd,
                  fill = "Ward"), alpha = 0.3) +
  geom_line(data = subset(data_regions_eng, region %in% c("South East", "South West", "East of England")),
            aes(date, Ward, colour = "Ward", linetype = "1Data"), size = 1) +
  geom_line(data = subset(res_regions_eng, region %in% c("South East", "South West", "East of England")),
            aes(date, CC, colour = "CC", linetype = "2Model"), size = 1) +
  geom_ribbon(data = subset(res_regions_eng, region %in% c("South East", "South West", "East of England")),
              aes(date,
                  ymin = CC - 1.96*CC_sd,
                  ymax = CC + 1.96*CC_sd,
                  fill = "CC"), alpha = 0.3) +
  geom_line(data = subset(data_regions_eng, region %in% c("South East", "South West", "East of England")),
            aes(date, CC, colour = "CC", linetype = "1Data"), size = 1) +
  scale_linetype_discrete(labels = c("Data", "Model")) +
  facet_wrap(~region) +
  guides(fill = F) +
  labs(x = "Date", y = "Bed occupancy", color = "Bed type:", linetype = "Value:") +
  theme_bw()

#all together

legend = get_legend(p1 + theme(legend.position = "bottom"))
plot_grid(p1 + theme(legend.position = "none"),
          p2 + theme(legend.position = "none"),
          p3 + theme(legend.position = "none"),
          legend,
          ncol = 1,
          rel_heights = c(1,1,1,0.2))


```
## By region (full pathways)

```{r region_plots_pathways, fig.width = 9, fig.height = 6, dpi = 300}

#london midlands 
p1 = ggplot() +
  geom_line(data = subset(res_regions_eng_pathways, region %in% c("London", "Midlands")),
            aes(date, Ward, colour = "Ward", linetype = "2Model"), size = 1) +
  geom_ribbon(data = subset(res_regions_eng_pathways, region %in% c("London", "Midlands")),
              aes(date,
                  ymin = Ward - 1.96*Ward_sd,
                  ymax = Ward + 1.96*Ward_sd,
                  fill = "Ward"), alpha = 0.3) +
  geom_line(data = subset(data_regions_eng, region %in% c("London", "Midlands")),
            aes(date, Ward, colour = "Ward", linetype = "1Data"), size = 1) +
  geom_line(data = subset(res_regions_eng_pathways, region %in% c("London", "Midlands")),
            aes(date, CC, colour = "CC", linetype = "2Model"), size = 1) +
  geom_ribbon(data = subset(res_regions_eng_pathways, region %in% c("London", "Midlands")),
              aes(date,
                  ymin = CC - 1.96*CC_sd,
                  ymax = CC + 1.96*CC_sd,
                  fill = "CC"), alpha = 0.3) +
  geom_line(data = subset(data_regions_eng, region %in% c("London", "Midlands")),
            aes(date, CC, colour = "CC", linetype = "1Data"), size = 1) +
  scale_linetype_discrete(labels = c("Data", "Model")) +
  facet_wrap(~region) +
  guides(fill = F) +
  labs(x = "Date", y = "Bed occupancy", color = "Bed type:", linetype = "Value:") +
  theme_bw()

#north east north west
p2 = ggplot() +
  geom_line(data = subset(res_regions_eng_pathways, region %in% c("North East and Yorkshire", "North West")),
            aes(date, Ward, colour = "Ward", linetype = "2Model"), size = 1) +
  geom_ribbon(data = subset(res_regions_eng_pathways, region %in% c("North East and Yorkshire", "North West")),
              aes(date,
                  ymin = Ward - 1.96*Ward_sd,
                  ymax = Ward + 1.96*Ward_sd,
                  fill = "Ward"), alpha = 0.3) +
  geom_line(data = subset(data_regions_eng, region %in% c("North East and Yorkshire", "North West")),
            aes(date, Ward, colour = "Ward", linetype = "1Data"), size = 1) +
  geom_line(data = subset(res_regions_eng_pathways, region %in% c("North East and Yorkshire", "North West")),
            aes(date, CC, colour = "CC", linetype = "2Model"), size = 1) +
  geom_ribbon(data = subset(res_regions_eng_pathways, region %in% c("North East and Yorkshire", "North West")),
              aes(date,
                  ymin = CC - 1.96*CC_sd,
                  ymax = CC + 1.96*CC_sd,
                  fill = "CC"), alpha = 0.3) +
  geom_line(data = subset(data_regions_eng, region %in% c("North East and Yorkshire", "North West")),
            aes(date, CC, colour = "CC", linetype = "1Data"), size = 1) +
  scale_linetype_discrete(labels = c("Data", "Model")) +
  facet_wrap(~region) +
  guides(fill = F) +
  labs(x = "Date", y = "Bed occupancy", color = "Bed type:", linetype = "Value:") +
  theme_bw()

#south east east south west
p3 = ggplot() +
  geom_line(data = subset(res_regions_eng_pathways, region %in% c("South East", "South West", "East of England")),
            aes(date, Ward, colour = "Ward", linetype = "2Model"), size = 1) +
  geom_ribbon(data = subset(res_regions_eng_pathways, region %in% c("South East", "South West", "East of England")),
              aes(date,
                  ymin = Ward - 1.96*Ward_sd,
                  ymax = Ward + 1.96*Ward_sd,
                  fill = "Ward"), alpha = 0.3) +
  geom_line(data = subset(data_regions_eng, region %in% c("South East", "South West", "East of England")),
            aes(date, Ward, colour = "Ward", linetype = "1Data"), size = 1) +
  geom_line(data = subset(res_regions_eng_pathways, region %in% c("South East", "South West", "East of England")),
            aes(date, CC, colour = "CC", linetype = "2Model"), size = 1) +
  geom_ribbon(data = subset(res_regions_eng_pathways, region %in% c("South East", "South West", "East of England")),
              aes(date,
                  ymin = CC - 1.96*CC_sd,
                  ymax = CC + 1.96*CC_sd,
                  fill = "CC"), alpha = 0.3) +
  geom_line(data = subset(data_regions_eng, region %in% c("South East", "South West", "East of England")),
            aes(date, CC, colour = "CC", linetype = "1Data"), size = 1) +
  scale_linetype_discrete(labels = c("Data", "Model")) +
  facet_wrap(~region) +
  guides(fill = F) +
  labs(x = "Date", y = "Bed occupancy", color = "Bed type:", linetype = "Value:") +
  theme_bw()

#all together

legend = get_legend(p1 + theme(legend.position = "bottom"))
plot_grid(p1 + theme(legend.position = "none"),
          p2 + theme(legend.position = "none"),
          p3 + theme(legend.position = "none"),
          legend,
          ncol = 1,
          rel_heights = c(1,1,1,0.2))


```



## Plot fit by region

```{r plot_fit_region, fig.width = 9, fig.height = 6, dpi = 300}

#london midlands 
p1 = ggplot() +
  geom_line(data = subset(fit_regions_eng, region %in% c("London", "Midlands")),
            aes(date, Ward, colour = "Ward", linetype = "2Model"), size = 1) +
  geom_ribbon(data = subset(fit_regions_eng, region %in% c("London", "Midlands")),
              aes(date,
                  ymin = Ward - 1.96*Ward_sd,
                  ymax = Ward + 1.96*Ward_sd,
                  fill = "Ward"), alpha = 0.3) +
  geom_line(data = subset(data_regions_eng, region %in% c("London", "Midlands")),
            aes(date, Ward, colour = "Ward", linetype = "1Data"), size = 1) +
  geom_line(data = subset(fit_regions_eng, region %in% c("London", "Midlands")),
            aes(date, CC, colour = "CC", linetype = "2Model"), size = 1) +
  geom_ribbon(data = subset(fit_regions_eng, region %in% c("London", "Midlands")),
              aes(date,
                  ymin = CC - 1.96*CC_sd,
                  ymax = CC + 1.96*CC_sd,
                  fill = "CC"), alpha = 0.3) +
  geom_line(data = subset(data_regions_eng, region %in% c("London", "Midlands")),
            aes(date, CC, colour = "CC", linetype = "1Data"), size = 1) +
  scale_linetype_discrete(labels = c("Data", "Model")) +
  facet_wrap(~region) +
  guides(fill = F) +
  labs(x = "Date", y = "Bed occupancy", color = "Bed type:", linetype = "Value:") +
  theme_bw()

#north east north west
p2 = ggplot() +
  geom_line(data = subset(fit_regions_eng, region %in% c("North East and Yorkshire", "North West")),
            aes(date, Ward, colour = "Ward", linetype = "2Model"), size = 1) +
  geom_ribbon(data = subset(fit_regions_eng, region %in% c("North East and Yorkshire", "North West")),
              aes(date,
                  ymin = Ward - 1.96*Ward_sd,
                  ymax = Ward + 1.96*Ward_sd,
                  fill = "Ward"), alpha = 0.3) +
  geom_line(data = subset(data_regions_eng, region %in% c("North East and Yorkshire", "North West")),
            aes(date, Ward, colour = "Ward", linetype = "1Data"), size = 1) +
  geom_line(data = subset(fit_regions_eng, region %in% c("North East and Yorkshire", "North West")),
            aes(date, CC, colour = "CC", linetype = "2Model"), size = 1) +
  geom_ribbon(data = subset(fit_regions_eng, region %in% c("North East and Yorkshire", "North West")),
              aes(date,
                  ymin = CC - 1.96*CC_sd,
                  ymax = CC + 1.96*CC_sd,
                  fill = "CC"), alpha = 0.3) +
  geom_line(data = subset(data_regions_eng, region %in% c("North East and Yorkshire", "North West")),
            aes(date, CC, colour = "CC", linetype = "1Data"), size = 1) +
  scale_linetype_discrete(labels = c("Data", "Model")) +
  facet_wrap(~region) +
  guides(fill = F) +
  labs(x = "Date", y = "Bed occupancy", color = "Bed type:", linetype = "Value:") +
  theme_bw()

#south east east south west
p3 = ggplot() +
  geom_line(data = subset(fit_regions_eng, region %in% c("South East", "South West", "East of England")),
            aes(date, Ward, colour = "Ward", linetype = "2Model"), size = 1) +
  geom_ribbon(data = subset(fit_regions_eng, region %in% c("South East", "South West", "East of England")),
              aes(date,
                  ymin = Ward - 1.96*Ward_sd,
                  ymax = Ward + 1.96*Ward_sd,
                  fill = "Ward"), alpha = 0.3) +
  geom_line(data = subset(data_regions_eng, region %in% c("South East", "South West", "East of England")),
            aes(date, Ward, colour = "Ward", linetype = "1Data"), size = 1) +
  geom_line(data = subset(fit_regions_eng, region %in% c("South East", "South West", "East of England")),
            aes(date, CC, colour = "CC", linetype = "2Model"), size = 1) +
  geom_ribbon(data = subset(fit_regions_eng, region %in% c("South East", "South West", "East of England")),
              aes(date,
                  ymin = CC - 1.96*CC_sd,
                  ymax = CC + 1.96*CC_sd,
                  fill = "CC"), alpha = 0.3) +
  geom_line(data = subset(data_regions_eng, region %in% c("South East", "South West", "East of England")),
            aes(date, CC, colour = "CC", linetype = "1Data"), size = 1) +
  scale_linetype_discrete(labels = c("Data", "Model")) +
  facet_wrap(~region) +
  guides(fill = F) +
  labs(x = "Date", y = "Bed occupancy", color = "Bed type:", linetype = "Value:") +
  theme_bw()

#all together

legend = get_legend(p1 + theme(legend.position = "bottom"))
plot_grid(p1 + theme(legend.position = "none"),
          p2 + theme(legend.position = "none"),
          p3 + theme(legend.position = "none"),
          legend,
          ncol = 1,
          rel_heights = c(1,1,1,0.2))


```


## Plot fit by region (proba)

```{r plot_fit_region_proba, fig.width = 9, fig.height = 6, dpi = 300}

#london midlands 
p1 = ggplot() +
  geom_line(data = subset(fit_regions_eng_proba, region %in% c("London", "Midlands")),
            aes(date, Ward, colour = "Ward", linetype = "2Model"), size = 1) +
  geom_ribbon(data = subset(fit_regions_eng_proba, region %in% c("London", "Midlands")),
              aes(date,
                  ymin = Ward - 1.96*Ward_sd,
                  ymax = Ward + 1.96*Ward_sd,
                  fill = "Ward"), alpha = 0.3) +
  geom_line(data = subset(data_regions_eng, region %in% c("London", "Midlands")),
            aes(date, Ward, colour = "Ward", linetype = "1Data"), size = 1) +
  geom_line(data = subset(fit_regions_eng_proba, region %in% c("London", "Midlands")),
            aes(date, CC, colour = "CC", linetype = "2Model"), size = 1) +
  geom_ribbon(data = subset(fit_regions_eng_proba, region %in% c("London", "Midlands")),
              aes(date,
                  ymin = CC - 1.96*CC_sd,
                  ymax = CC + 1.96*CC_sd,
                  fill = "CC"), alpha = 0.3) +
  geom_line(data = subset(data_regions_eng, region %in% c("London", "Midlands")),
            aes(date, CC, colour = "CC", linetype = "1Data"), size = 1) +
  scale_linetype_discrete(labels = c("Data", "Model")) +
  facet_wrap(~region) +
  guides(fill = F) +
  labs(x = "Date", y = "Bed occupancy", color = "Bed type:", linetype = "Value:") +
  theme_bw()

#north east north west
p2 = ggplot() +
  geom_line(data = subset(fit_regions_eng_proba, region %in% c("North East and Yorkshire", "North West")),
            aes(date, Ward, colour = "Ward", linetype = "2Model"), size = 1) +
  geom_ribbon(data = subset(fit_regions_eng_proba, region %in% c("North East and Yorkshire", "North West")),
              aes(date,
                  ymin = Ward - 1.96*Ward_sd,
                  ymax = Ward + 1.96*Ward_sd,
                  fill = "Ward"), alpha = 0.3) +
  geom_line(data = subset(data_regions_eng, region %in% c("North East and Yorkshire", "North West")),
            aes(date, Ward, colour = "Ward", linetype = "1Data"), size = 1) +
  geom_line(data = subset(fit_regions_eng_proba, region %in% c("North East and Yorkshire", "North West")),
            aes(date, CC, colour = "CC", linetype = "2Model"), size = 1) +
  geom_ribbon(data = subset(fit_regions_eng_proba, region %in% c("North East and Yorkshire", "North West")),
              aes(date,
                  ymin = CC - 1.96*CC_sd,
                  ymax = CC + 1.96*CC_sd,
                  fill = "CC"), alpha = 0.3) +
  geom_line(data = subset(data_regions_eng, region %in% c("North East and Yorkshire", "North West")),
            aes(date, CC, colour = "CC", linetype = "1Data"), size = 1) +
  scale_linetype_discrete(labels = c("Data", "Model")) +
  facet_wrap(~region) +
  guides(fill = F) +
  labs(x = "Date", y = "Bed occupancy", color = "Bed type:", linetype = "Value:") +
  theme_bw()

#south east east south west
p3 = ggplot() +
  geom_line(data = subset(fit_regions_eng_proba, region %in% c("South East", "South West", "East of England")),
            aes(date, Ward, colour = "Ward", linetype = "2Model"), size = 1) +
  geom_ribbon(data = subset(fit_regions_eng_proba, region %in% c("South East", "South West", "East of England")),
              aes(date,
                  ymin = Ward - 1.96*Ward_sd,
                  ymax = Ward + 1.96*Ward_sd,
                  fill = "Ward"), alpha = 0.3) +
  geom_line(data = subset(data_regions_eng, region %in% c("South East", "South West", "East of England")),
            aes(date, Ward, colour = "Ward", linetype = "1Data"), size = 1) +
  geom_line(data = subset(fit_regions_eng_proba, region %in% c("South East", "South West", "East of England")),
            aes(date, CC, colour = "CC", linetype = "2Model"), size = 1) +
  geom_ribbon(data = subset(fit_regions_eng_proba, region %in% c("South East", "South West", "East of England")),
              aes(date,
                  ymin = CC - 1.96*CC_sd,
                  ymax = CC + 1.96*CC_sd,
                  fill = "CC"), alpha = 0.3) +
  geom_line(data = subset(data_regions_eng, region %in% c("South East", "South West", "East of England")),
            aes(date, CC, colour = "CC", linetype = "1Data"), size = 1) +
  scale_linetype_discrete(labels = c("Data", "Model")) +
  facet_wrap(~region) +
  guides(fill = F) +
  labs(x = "Date", y = "Bed occupancy", color = "Bed type:", linetype = "Value:") +
  theme_bw()

#all together

legend = get_legend(p1 + theme(legend.position = "bottom"))
plot_grid(p1 + theme(legend.position = "none"),
          p2 + theme(legend.position = "none"),
          p3 + theme(legend.position = "none"),
          legend,
          ncol = 1,
          rel_heights = c(1,1,1,0.2))


```

## Massive figure for both datasets, LoS flipped between datasets

```{r eng_res_melted_uclh}

eng_hospital_data = hospital_data %>%
  filter(geography == "England")

Mres_average_LoS = res_average_LoS_eng_uclh %>%
  select(-time) %>%
  melt(id.vars = "date")
Mres_average_LoS$grouping = "Model - average LoS"

Mres_full_pathways = res_full_pathways_eng_uclh %>%
  select(-time) %>%
  melt(id.vars = "date")
Mres_full_pathways$grouping = "Model - full pathways"

Mres_all = rbind(Mres_average_LoS, Mres_full_pathways)

#separate mean and sd
results_mean = Mres_all %>%
  filter(!grepl("sd", variable)) %>%
  arrange(variable) %>%
  mutate(date = as.Date(date))

results_sd = Mres_all %>% 
  filter(grepl("sd", variable)) %>%
  arrange(variable) %>%
  mutate(date = as.Date(date))

results_sd$min = results_mean$value - 1.96*results_sd$value
results_sd$max = results_mean$value + 1.96*results_sd$value


```

```{r uclh_res_melted_cocin}

Mres_average_LoS_uclh = res_average_LoS_uclh_cocin %>%
  select(-time) %>%
  melt(id.vars = "date")
Mres_average_LoS_uclh$grouping = "Model - average LoS"

Mres_full_pathways_uclh = res_full_pathways_uclh_cocin %>%
  select(-time) %>%
  melt(id.vars = "date")
Mres_full_pathways_uclh$grouping = "Model - full pathways"

Mres_all_uclh = rbind(Mres_average_LoS_uclh, Mres_full_pathways_uclh)

#separate mean and sd
results_mean_uclh = Mres_all_uclh %>%
  filter(!grepl("sd", variable)) %>%
  arrange(variable) %>%
  mutate(date = as.Date(date))

results_sd_uclh = Mres_all_uclh %>% 
  filter(grepl("sd", variable)) %>%
  arrange(variable) %>%
  mutate(date = as.Date(date))

results_sd_uclh$min = results_mean_uclh$value - 1.96*results_sd_uclh$value
results_sd_uclh$max = results_mean_uclh$value + 1.96*results_sd_uclh$value


```


```{r massiveplot_flippedLoS, fig.width = 12, fig.height = 7, dpi = 300}

#create plot A
plotA = ggplot() +
  geom_line(data = results_mean_uclh %>% filter(grouping == "Model - average LoS"),
            aes(date, value, colour = variable, linetype = "Model"), size = 1) +
  geom_ribbon(data = subset(results_sd_uclh, variable == "CC_sd") %>%
                filter(grouping == "Model - average LoS"), 
              aes(date,
                  ymin = min,
                  ymax = max,
                  fill = "CC"),
              alpha = 0.3) +
  geom_ribbon(data = subset(results_sd_uclh, variable == "Ward_sd")%>%
                filter(grouping == "Model - average LoS"),
              aes(date,
                  ymin = min,
                  ymax = max,
                  fill = "Ward"),
              alpha = 0.3) +
  geom_line(data = uclh_data %>% filter(value_type == "Ward"),
            aes(date, value, linetype = "Data", colour = "Ward"), size = 1) +
  geom_line(data = uclh_data %>% filter(value_type == "CC"),
            aes(date, value, linetype = "Data", colour = "CC"), size = 1) +
  labs(colour = "Bed type:", linetype = "Value:",
       x = "Time (days)", y = "Bed occupancy") +
  guides(fill = FALSE) +
  scale_x_date(breaks = "1 week", date_labels = format("%d %b %Y")) +
  theme_bw()

#create plot B
plotB = ggplot() +
  geom_line(data = results_mean_uclh %>% filter(grouping == "Model - full pathways"),
            aes(date, value, colour = variable, linetype = "Model"), size = 1) +
  geom_ribbon(data = subset(results_sd_uclh, variable == "CC_sd") %>%
                filter(grouping == "Model - full pathways"), 
              aes(date,
                  ymin = min,
                  ymax = max,
                  fill = "CC"),
              alpha = 0.3) +
  geom_ribbon(data = subset(results_sd_uclh, variable == "Ward_sd")%>%
                filter(grouping == "Model - full pathways"),
              aes(date,
                  ymin = min,
                  ymax = max,
                  fill = "Ward"),
              alpha = 0.3) +
  geom_line(data = uclh_data %>% filter(value_type == "Ward"),
            aes(date, value, linetype = "Data", colour = "Ward"), size = 1) +
  geom_line(data = uclh_data %>% filter(value_type == "CC"),
            aes(date, value, linetype = "Data", colour = "CC"), size = 1) +
  labs(colour = "Bed type:", linetype = "Value:",
       x = "Time (days)", y = "Bed occupancy") +
  guides(fill = FALSE) +
  scale_x_date(breaks = "1 week", date_labels = format("%d %b %Y")) +
  theme_bw()

#create plot C
plotC = ggplot() +
  geom_line(data = results_mean %>% filter(grouping == "Model - average LoS"),
            aes(date, value, colour = variable, linetype = "Model"), size = 1) +
  geom_ribbon(data = subset(results_sd, variable == "CC_sd") %>%
                filter(grouping == "Model - average LoS"), 
              aes(date,
                  ymin = min,
                  ymax = max,
                  fill = "CC"),
              alpha = 0.3) +
  geom_ribbon(data = subset(results_sd, variable == "Ward_sd")%>%
                filter(grouping == "Model - average LoS"),
              aes(date,
                  ymin = min,
                  ymax = max,
                  fill = "Ward"),
              alpha = 0.3) +
  geom_line(data = eng_hospital_data %>% filter(value_type == "hospital_prev"),
            aes(date, value, linetype = "Data", colour = "Ward"), size = 1) +
  geom_line(data = eng_hospital_data %>% filter(value_type == "icu_prev"),
            aes(date, value, linetype = "Data", colour = "CC"), size = 1) +
  labs(colour = "Bed type:", linetype = "Value:",
       x = "Time (days)", y = "Bed occupancy") +
  guides(fill = FALSE) +
  scale_x_date(breaks = "1 month", date_labels = format("%d %b %Y")) +
  theme_bw()

#create plot D
plotD = ggplot() +
  geom_line(data = results_mean %>% filter(grouping == "Model - full pathways"),
            aes(date, value, colour = variable, linetype = "Model"), size = 1) +
  geom_ribbon(data = subset(results_sd, variable == "CC_sd") %>%
                filter(grouping == "Model - full pathways"), 
              aes(date,
                  ymin = min,
                  ymax = max,
                  fill = "CC"),
              alpha = 0.3) +
  geom_ribbon(data = subset(results_sd, variable == "Ward_sd")%>%
                filter(grouping == "Model - full pathways"),
              aes(date,
                  ymin = min,
                  ymax = max,
                  fill = "Ward"),
              alpha = 0.3) +
  geom_line(data = eng_hospital_data %>% filter(value_type == "hospital_prev"),
            aes(date, value, linetype = "Data", colour = "Ward"), size = 1) +
  geom_line(data = eng_hospital_data %>% filter(value_type == "icu_prev"),
            aes(date, value, linetype = "Data", colour = "CC"), size = 1) +
  labs(colour = "Bed type:", linetype = "Value:",
       x = "Time (days)", y = "Bed occupancy") +
  guides(fill = FALSE) +
  scale_x_date(breaks = "1 month", date_labels = format("%d %b %Y")) +
  theme_bw()


legend = get_legend(plotA + theme(legend.position = "bottom"))
plot_grid(
  plot_grid(
    plotA+theme(legend.position = "none"), plotB+theme(legend.position = "none"),
    plotC+theme(legend.position = "none"), plotD+theme(legend.position = "none"),
    labels = c("A)", "B)", "C)", "D)"), 
    label_size = 12, hjust = -0.3
  ),
  legend,
  nrow = 2,
  rel_heights = c(1, 0.1))
```


## Supplementary table - sq diff 


```{r sqdiffbyreg}

sqdiff_regions = data.frame(region = unique(data_regions_eng$region),
                            sqdiff_average = 0,
                            sqdiff_pathways = 0,
                            sqdiff_los_fit = 0,
                            sqdiff_cc_fit = 0)

for(i in 1:nrow(sqdiff_regions)){
  
  region_i = sqdiff_regions$region[i]
  
  dat_i = data_regions_eng %>%
    filter(region == region_i)
  
  av_i = res_regions_eng %>%
    filter(region == region_i)
  path_i = res_regions_eng_pathways %>%
    filter(region == region_i)
  los_fit_i = fit_regions_eng %>%
    filter(region == region_i)
  cc_fit_i = fit_regions_eng_proba %>%
    filter(region == region_i)
  
  sqdiff_regions$sqdiff_average[i] = round(sum((dat_i$Ward-av_i$Ward)^2,
                                               (dat_i$CC-av_i$CC)^2, na.rm = T)/10^5, 2)
  sqdiff_regions$sqdiff_pathways[i] = round(sum((dat_i$Ward-path_i$Ward)^2,
                                                (dat_i$CC-path_i$CC)^2, na.rm = T)/10^5, 2)
  sqdiff_regions$sqdiff_los_fit[i] = round(sum((dat_i$Ward-los_fit_i$Ward)^2,
                                               (dat_i$CC-los_fit_i$CC)^2, na.rm = T)/10^5, 2)
  sqdiff_regions$sqdiff_cc_fit[i] = round(sum((dat_i$Ward-cc_fit_i$Ward)^2,
                                              (dat_i$CC-cc_fit_i$CC)^2, na.rm = T)/10^5, 2)
  
}


kable(sqdiff_regions, format = "latex", booktabs = T,
      caption = "Squared difference between model predicted occupancy and data. CC: critical care (ICU or HDU).",
      col.names = c("NHS region", "Average LoS",
                    "Full pathways", "Best fit LoS", "Best fit prop to CC")) %>%
  kable_styling(latex_options = c("striped"))

```