WPI_aggregate_desc.qmd

---
title: "NHS England Community stategy"
author: "Alexander Lawless"
date: last-modified
date-format: "DD-MM-YYYY"
title-block-banner: "#f9bf07"
title-block-banner-color: "#333739"
format:
  html:
    embed-resources: true
    smooth-scroll: true
    theme: cosmo
    fontcolor: black
    toc: true
    toc-location: left
    toc-title: Summary
    toc-depth: 3
editor: visual
css: styles.css
---

```{r}
#| label: setup
#| include: false
#| echo: false
#| warning: false

library(tidyverse)
library(janitor)
library(DT)
library(sf)
library(readxl)
library(patchwork)
library(GGally)
library(plotly)
library(reshape2)
library(ggrepel)

knitr::opts_chunk$set(echo = FALSE, warning=FALSE, message=FALSE, fig.width =12, fig.height = 9)

setwd("C:/Users/alexander.lawless/OneDrive - Midlands and Lancashire CSU/Work/1. Projects/2024_25/Community services modelling/r_repos/Community_Strategies")

# Functions/lookups ----
create_dt <- function(x) {
    
    DT::datatable(
        x
        , extensions = "Buttons"
        , options = list(
            dom = "Blfrtip"
            , buttons = c("copy", "csv")
            , lengthMenu = list(
                c(10, 25, 50, -1)
                , c(10, 25, 50, "All")
                ),
            width = "400px",  # Specify the width
            height = "200px" # Specify the height
            )
        )
}

ethnicity_lookup <- tribble(
  ~Code, ~Description,
  "A", "White - British",
  "B", "White - Irish",
  "C", "White - Any other White background",
  "D", "Mixed - White and Black Caribbean",
  "E", "Mixed - White and Black African",
  "F", "Mixed - White and Asian",
  "G", "Mixed - Any other mixed background",
  "H", "Asian or Asian British - Indian",
  "J", "Asian or Asian British - Pakistani",
  "K", "Asian or Asian British - Bangladeshi",
  "L", "Asian or Asian British - Any other Asian background",
  "M", "Black or Black British - Caribbean",
  "N", "Black or Black British - African",
  "P", "Black or Black British - Any other Black background",
  "R", "Other Ethnic Groups - Chinese",
  "S", "Other Ethnic Groups - Any other ethnic group",
  "Z", "Not stated"
)

# Set SU theme ####
SU_colours <- c(
  `orange`                     = grDevices::rgb(248, 191, 7, maxColorValue = 255), # "#f9bf07",
  `charcoal`                   = grDevices::rgb(44, 40, 37, maxColorValue = 255), # "#2c2825",
  `slate`                      = grDevices::rgb(104, 111, 115, maxColorValue = 255), # "#686f73",
  `blue`                       = grDevices::rgb(88, 29, 193, maxColorValue = 255), # "#5881c1",
  `red`                        = grDevices::rgb(236, 101, 85, maxColorValue = 255), # "#ec6555",
  # additional accent colours from word doc template
  `yellow`                     = grDevices::rgb(252, 229, 155, maxColorValue = 255),
  `grey`                       = grDevices::rgb(163, 168, 172, maxColorValue = 255),
  `white`                      = grDevices::rgb(255, 255, 255, maxColorValue = 255),
  # light and dark ends from colour theme in word doc
  `light orange`               = grDevices::rgb(253, 242, 205, maxColorValue = 255),
  `dark orange`                = grDevices::rgb(124, 95, 3, maxColorValue = 255),
  `light charcoal`             = grDevices::rgb(235, 233, 231, maxColorValue = 255),
  `dark charcoal`              =        "#000000", # black
  `light slate`                = grDevices::rgb(224, 226, 227, maxColorValue = 255),
  `dark slate`                 = grDevices::rgb(51, 55, 57, maxColorValue = 255),
  `light blue`                 = grDevices::rgb(221, 229, 242, maxColorValue = 255),
  `dark blue`                  = grDevices::rgb(38, 61, 102, maxColorValue = 255),
  `light red`                  = grDevices::rgb(251, 224, 220, maxColorValue = 255),
  `dark red`                   = grDevices::rgb(144, 29, 16, maxColorValue = 255),
  `light yellow`               = grDevices::rgb(254, 249, 235, maxColorValue = 255),
  `dark yellow`                = grDevices::rgb(197, 152, 5, maxColorValue = 255),
  `light grey`                 = grDevices::rgb(236, 237, 238, maxColorValue = 255),
  `dark grey`                  = grDevices::rgb(79, 84, 88, maxColorValue = 255),
  `light white`                = grDevices::rgb(242, 242, 242, maxColorValue = 255),
  `dark white`                 = grDevices::rgb(127, 127, 127, maxColorValue = 255),
  `red2`                       = grDevices::rgb(215, 25, 28, maxColorValue = 255),
  `orange2`                    = grDevices::rgb(253, 174, 97, maxColorValue = 255),
  `yellow2`                    = grDevices::rgb(255, 255, 191, maxColorValue = 255),
  `green2`                     = grDevices::rgb(171, 221, 164, maxColorValue = 255),
  `blue2`                      = grDevices::rgb(43, 131, 186, maxColorValue = 255) # "#2b83ba"
)

SU_cols <- function(...) {
  cols <- c(...)

  if (is.null(cols)) {
    return(SU_colours)
  }

  SU_colours[cols]
}

SU_palettes <- list(
  `main` = SU_cols("orange", "charcoal", "slate", "blue", "red"),
  `oranges` = SU_cols("light orange", "orange", "dark orange"),
  `slates` = SU_cols("light slate", "slate", "dark slate"),
  `mixed` = SU_cols("dark red", "orange", "yellow", "light blue", "slate"),
  `oj_coal` = SU_cols("yellow", "orange", "red", "dark red", "dark charcoal"),
  `oj_red` = SU_cols("yellow", "orange", "red", "dark red"),
  `white_oj_coal` = SU_cols("white", "yellow", "orange", "red", "dark red", "dark charcoal"), # added since shared
  `lyellow_oj_coal` = SU_cols("light yellow", "orange", "red", "dark red", "dark charcoal"), # added since shared
  `wy_oj_coal` = SU_cols("white", "light yellow", "yellow", "orange", "red", "dark red", "charcoal", "dark charcoal"),
  `red_coal` = SU_cols("red", "dark red", "charcoal", "dark charcoal"),
  `blue_yellow_red` = SU_cols("red2", "orange2", "yellow2", "green2", "blue2"),
  `red_yellow_blue` = SU_cols("blue2", "green2", "yellow2", "orange2", "red2")
)


SU_pal <- function(palette = "main", reverse = FALSE, ...) {
  pal <- SU_palettes[[palette]]

  if (reverse) pal <- rev(pal)

  colorRampPalette(pal, ...)
}


scale_color_SU <- function(palette = "main", discrete = TRUE, reverse = FALSE, ...) {
  pal <- SU_pal(palette = palette, reverse = reverse)

  if (discrete) {
    discrete_scale("colour", paste0("SU_", palette), palette = pal, ...)
  } else {
    scale_color_gradientn(colours = pal(256), ...)
  }
}

scale_fill_SU <- function(palette = "main", discrete = TRUE, reverse = FALSE, ...) {
  pal <- SU_pal(palette = palette, reverse = reverse)

  if (discrete) {
    discrete_scale("fill", paste0("SU_", palette), palette = pal, ...)
  } else {
    scale_fill_gradientn(colours = pal(256), ...)
  }
}

theme_SU <- function(base_size) {
  theme_minimal(
    # base_family = "Segoe UI",
    base_size = 12
  ) %+replace%
    theme(
      axis.title = element_text(size = 11, face = "bold", colour = SU_cols("charcoal")),
      plot.title = element_text(hjust = 0, face = "bold", size = 12, colour = SU_cols("charcoal"), margin = margin(b = 4, unit = "pt")),
      plot.subtitle = element_text(hjust = 0, face = "italic", size = 10, colour = SU_cols("charcoal"), margin = margin(b = 4, unit = "pt")),
      plot.caption = element_text(hjust = 0, face = "italic", size = 9, colour = SU_cols("slate"), margin = margin(b = 4, unit = "pt")),
      legend.text = element_text(size = 10, colour = SU_cols("charcoal")),
      legend.title = element_text(face = "bold", size = 11, colour = SU_cols("charcoal"), margin = margin(b = 4, unit = "pt"))
    )
}

theme_set(theme_SU())

# Access comm strat grouped data ----

aggregate_data <- 
  read_csv("aggregate_data_recombined.csv") |> 
  mutate(flag_der_falls = 
           case_when(!is.na(flag_fall_imp_frac) ~ flag_fall_imp_frac,
                     TRUE ~ flag_fall_imp_tend)) |> 
  filter(der_activity_month <= 202408)

aggregate_data_core_cohorts <-
  aggregate_data |> 
  filter(!(is.na(flag_frail) & is.na(flag_eol) & is.na(flag_elderly_emergency) & is.na(flag_der_falls)))


# Split data into groups of interest ----

aggregate_data_frail <-
  aggregate_data %>% 
  drop_na(flag_frail) %>% 
  mutate(year  = lubridate::year(month))

aggregate_data_eol <-
  aggregate_data %>% 
  drop_na(flag_eol) %>% 
  mutate(year  = lubridate::year(month))

aggregate_data_falls <-
  aggregate_data %>% 
  drop_na(flag_der_falls) %>% 
  mutate(year  = lubridate::year(month))

#aggregate_data_falls <-
#  aggregate_data %>% 
#  filter(!is.na(flag_falls_exp) |
#           !is.na(flag_fall_imp_frac) | 
#           !is.na(flag_fall_imp_tend)) %>% 
#  mutate(year  = lubridate::year(month))

aggregate_data_elderly_emergency <-
  aggregate_data %>% 
  drop_na(flag_elderly_emergency) %>% 
  mutate(year  = lubridate::year(month))

aggregate_data_amb_chronic <-
  aggregate_data %>% 
  drop_na(amb_chronic) %>% 
  mutate(year  = lubridate::year(month))

aggregate_data_amb_acute <-
  aggregate_data %>% 
  drop_na(amb_acute) %>% 
  mutate(year  = lubridate::year(month))

aggregate_data_amb_vacc_prev <-
  aggregate_data %>% 
  drop_na(amb_vacc_prev) %>% 
  mutate(year  = lubridate::year(month))

aggregate_data_eol_1_year <-
  aggregate_data %>% 
  drop_na(death_location_type) %>% 
  mutate(year  = lubridate::year(month))

# Mapping ---- 

# Read in icb lookup
icb_lookup <- 
  read_csv("mapping/S_ICB_ICB_to_NHSE_R.csv") %>% 
  clean_names() %>% 
  select(icb23cd, icb23cdh, icb23nm, nhser23nm) %>% 
  distinct()

# Import icb shp file from geoportal
icb_23_shp <- 
  st_read("https://services1.arcgis.com/ESMARspQHYMw9BZ9/arcgis/rest/services/Integrated_Care_Boards_April_2023_EN_BGC/FeatureServer/0/query?outFields=*&where=1%3D1&f=geojson") %>% 
  clean_names()

# Age and sex adjusted population rate ----

icb_pop_2023 <- 
  read_excel("pop_estimates/sapehealthgeogstablefinal.xlsx", 
             sheet = "Mid-2022 ICB 2023", skip = 3) %>% 
  clean_names()

# Clean and group by ICB, sex and age-range
icb_pop_2023_sex_age_range <- 
  icb_pop_2023 %>%
  pivot_longer(cols = starts_with("f") | starts_with("m"), 
               names_to = "age_sex", 
               values_to = "population") %>%
  mutate(sex = as.numeric(ifelse(grepl("^m", age_sex), "1", "2")),
         age = as.numeric(sub("^[mf]", "", age_sex))) %>%
  mutate(age_range = 
           case_when(
             age >= 0 & age <= 4 ~ "0-4",
             age >= 5 & age <= 9 ~ "5-9",
             age >= 10 & age <= 14 ~ "10-14",
             age >= 15 & age <= 19 ~ "15-19",
             age >= 20 & age <= 24 ~ "20-24",
             age >= 25 & age <= 29 ~ "25-29",
             age >= 30 & age <= 34 ~ "30-34",
             age >= 35 & age <= 39 ~ "35-39",
             age >= 40 & age <= 44 ~ "40-44",
             age >= 45 & age <= 49 ~ "45-49",
             age >= 50 & age <= 54 ~ "50-54",
             age >= 55 & age <= 59 ~ "55-59",
             age >= 60 & age <= 64 ~ "60-64",
             age >= 65 & age <= 69 ~ "65-69",
             age >= 70 & age <= 74 ~ "70-74",
             age >= 75 & age <= 79 ~ "75-79",
             age >= 80 & age <= 84 ~ "80-84",
             age >= 85 & age <= 89 ~ "85-89",
             age >= 90 & age <= 94 ~ "90-94",
             age >= 95 & age <= 99 ~ "95-99",
             age >= 100 & age <= 104 ~ "100-104",
             age >= 105 ~ "105+",
             TRUE ~ NA
             )) %>% 
  #mutate(age_range = 
  #         cut(age, 
  #             breaks = seq(0, 90, by = 5), 
  #             right = FALSE, 
  #             labels = paste(seq(0, 85, by = 5))
  #             )
  #       ) %>% 
  group_by(icb_2023_code, icb_2023_name, sex, age_range) %>% 
  summarise(population = sum(population, na.rm= TRUE)) %>%
  ungroup() %>% 
  mutate(icb_2023_name = str_replace_all(icb_2023_name, "Integrated Care Board", "ICB")) %>% 
  mutate(icb_name_short_trans = 
           case_when(
             icb_2023_name == "NHS Bath and North East Somerset, Swindon and Wiltshire ICB" ~ "NHS Bath & NE S'set, S'don & W ICB",
             icb_2023_name == "NHS Bedfordshire, Luton and Milton Keynes ICB" ~ "NHS Beds, Luton & Milton Keynes ICB",
             icb_2023_name == "NHS Birmingham and Solihull ICB" ~ "NHS Birmingham And Solihull ICB",
             icb_2023_name == "NHS Bristol, North Somerset and South Gloucestershire ICB" ~ "NHS Bristol, N S'set & S Gloucs ICB",
             icb_2023_name == "NHS Buckinghamshire, Oxfordshire and Berkshire West ICB" ~ "NHS Bucks, Oxford & Berkshire W ICB",
             icb_2023_name == "NHS Cambridgeshire and Peterborough ICB" ~ "NHS Cambs & Peterborough ICB",
             icb_2023_name == "NHS Cornwall and the Isles of Scilly ICB" ~ "NHS Cornwall & Isles Of Scilly ICB",
             icb_2023_name == "NHS Herefordshire and Worcestershire ICB" ~ "NHS Herefords & Worcestershire ICB",
             icb_2023_name == "NHS Hertfordshire and West Essex ICB" ~ "NHS Hertfordshire & West Essex ICB",
             icb_2023_name == "NHS Lancashire and South Cumbria ICB" ~ "NHS Lancashire & South Cumbria ICB",
             icb_2023_name == "NHS Leicester, Leicestershire and Rutland ICB" ~ "NHS Leics, Leic'shire & Rutland ICB",
             icb_2023_name == "NHS North East and North Cumbria ICB" ~ "NHS North East & North Cumbria ICB",
             icb_2023_name == "NHS Nottingham and Nottinghamshire ICB" ~ "NHS Nott'ham & Nottinghamshire ICB",
             icb_2023_name == "NHS Shropshire, Telford and Wrekin ICB" ~ "NHS Shrops, Telford & Wrekin ICB",
             icb_2023_name == "NHS Staffordshire and Stoke-on-Trent ICB" ~ "NHS Staffordshire & Stoke-On-Tr ICB",
             icb_2023_name == "NHS Suffolk and North East Essex ICB" ~ "NHS Suffolk & North East Essex ICB",
             TRUE ~ icb_2023_name
           ))


# Apply standardisation to cohort sub-groups

standardisation_function <- function(data_input) {
  
crude_spell_rate <-
  data_input %>%  
  mutate(year = year(month),
         age_range = as.character(age_range)) %>% 
  group_by(year, icb_name_short, sex, age_range) %>% 
  summarise(spells = sum(spells)) %>%
  mutate(icb_code = str_extract(icb_name_short, "^[^:]+")) %>% 
  left_join(icb_lookup, by = c("icb_code" = "icb23cdh")) %>% 
  left_join(icb_pop_2023_sex_age_range %>% 
              mutate(age_range = as.character(age_range)),
            by = c("icb23cd" = "icb_2023_code", 
                   "sex", "age_range"), keep = F) %>%
  mutate(spell_rate = spells/population)

# Define a standard population (e.g., the total population in the dataset)
standard_population <- 
  icb_pop_2023_sex_age_range %>%
  group_by(age_range, sex) %>%
  summarise(population = sum(population))

# Merge with the standard population and calculate adjusted rate
adjusted_rate_icb <-
  crude_spell_rate %>%
  left_join(standard_population, by = c("age_range", "sex"), suffix = c("", "_standard")) %>% 
  drop_na(spell_rate) %>% 
  group_by(icb23cd, icb_name_short, year) %>%
  summarise(adjusted_spell_rate = sum(spell_rate * population_standard) / sum(population_standard) * 1000)

adjusted_rate_icb

}

adjusted_rate_sub_cohorts <-
  standardisation_function(aggregate_data) %>% 
  mutate(id = "rate") %>% 
  union_all(standardisation_function(aggregate_data_frail) %>% 
              mutate(id = "frail_rate")) %>% 
  union_all(standardisation_function(aggregate_data_falls) %>% 
              mutate(id = "falls_rate")) %>%
  union_all(standardisation_function(aggregate_data_eol) %>% 
              mutate(id = "eol_rate")) %>% 
  union_all(standardisation_function(aggregate_data_elderly_emergency) %>% 
              mutate(id = "elderly_emergency_rate")) %>% 
  ungroup() %>% 
  mutate(id_clean = case_when(id == "frail_rate" ~ "A. Frail",
                              id == "falls_rate" ~ "C. Falls",
                              id == "eol_rate" ~ "D. End of life",
                              id == "elderly_emergency_rate" ~ "B. Emergency elderly"))

crude_spell_rate <-
  aggregate_data %>%  
  mutate(year = year(month),
         age_range = as.character(age_range)) %>% 
  group_by(year, icb_name_short, sex, age_range) %>% 
  summarise(spells = sum(spells)) %>%
  mutate(icb_code = str_extract(icb_name_short, "^[^:]+")) %>% 
  left_join(icb_lookup, by = c("icb_code" = "icb23cdh")) %>% 
  left_join(icb_pop_2023_sex_age_range %>% 
              mutate(age_range = as.character(age_range)),
            by = c("icb23cd" = "icb_2023_code", 
                   "sex", "age_range"), keep = F) %>%
  mutate(spell_rate = spells/population)

# Avg LoS by ICB
avg_los_icb_cohort <- 
  read_csv("avg_los_icb_cohort.csv") %>% 
  clean_names()

# Create pairs plot
compare_admission_rates <-
  adjusted_rate_sub_cohorts |> 
  filter(year == 2023) |> 
  select(-icb23cd, -year, -id) |> 
  pivot_wider(id_cols = icb_name_short,
              names_from = id_clean,
              values_from = adjusted_spell_rate) |> 
  select(icb_name_short,`A. Frail`, `B. Emergency elderly`, `C. Falls`, `D. End of life`) 

thresholds <- sapply(compare_admission_rates[, 2:5], function(x) quantile(x, 0.9, na.rm = T))


# Admission rate map 
plot_function <- function(cohort, label) {
  
  icb_23_shp %>%
    select(icb23cd) %>%
    left_join(adjusted_rate_sub_cohorts %>%  
                filter(id == cohort,
                       year == 2023) %>% 
                mutate(quintile = cut(adjusted_spell_rate, 
                                      breaks = quantile(adjusted_spell_rate, probs = seq(0, 1, 0.2), na.rm = TRUE), 
                                      include.lowest = TRUE, 
                                      labels = FALSE)),
              by = "icb23cd") %>%
    
    ggplot() +
    geom_sf(aes(fill = factor(quintile))) +
    scale_fill_brewer(palette = "YlOrRd", 
                      name = paste0(label, ": Quintiles"),
                      labels = c("1st", "2nd", "3rd", "4th", "5th")) +
    theme(panel.grid = element_blank(),
          axis.text = element_blank()) 
}

frail_map <- plot_function("frail_rate", "A. Frail")
falls_map <- plot_function("falls_rate", "C. Falls")
eol_map <- plot_function("eol_rate", "D. EoL")
elderly_emergency_map <- plot_function("elderly_emergency_rate", "B. Elderly emergency")

# Plot patchwork map
combined_map_icb_rate_cohorts <-
  frail_map +
  elderly_emergency_map +
  falls_map +
  eol_map +
  
  plot_annotation(
    title = 'Mitigable activity by ICB - admission rate',
    subtitle = 'SUS 2023 | Admission identified: Frail, emergency elderly, falls & end of life',
    caption = 'Note: 1st quintile represents the lowest 20% admission rates'
  )

# Create function to draw funell plots
funnel_plot_function <- function(cohort, subtitle_text) {
  
  data <-
    adjusted_rate_sub_cohorts |> 
    filter(id_clean == {{cohort}}, 
           year == 2023) |> 
    left_join(icb_pop_2023_sex_age_range |> 
                group_by(icb_2023_code) |> 
                summarise(population = sum(population)),
              by = c("icb23cd" = "icb_2023_code")
              ) 
  
  mean_rate <- mean(data$adjusted_spell_rate)
  sd_rate <- sd(data$adjusted_spell_rate)
  
  # Create the funnel plot
  data |>
    mutate(std_from_mean = (adjusted_spell_rate - mean_rate) / sd_rate) |> 
    mutate(fill_text = 
             case_when(std_from_mean >= 1 ~ "1. Above 1 SD",
                       std_from_mean <= -1 ~ "3. Below 1 SD",
                       TRUE ~ "2. Within 1 SD")) |>
    mutate(icb_name_clean = str_sub(icb_name_short, 10,100)) |> 
    mutate(icb_name_clean = str_remove_all(icb_name_clean, " ICB"))|> 
    mutate(label = case_when(fill_text != "2. Within 1 SD" ~ icb_name_clean)) |> 
    
    ggplot(aes(x = population, y = adjusted_spell_rate, 
               colour = fill_text, 
               alpha = fill_text
    )) +
    geom_point(size = 4) +
    geom_label_repel(aes(label = label), size = 5, show.legend = FALSE) +
    geom_hline(yintercept = mean_rate, color = "blue", linetype = "dashed", linewidth = 1) +
    geom_hline(yintercept = mean_rate + sd_rate, color = "red", linetype = "dotted", linewidth = 1) +
    geom_hline(yintercept = mean_rate - sd_rate, color = "red", linetype = "dotted", linewidth = 1) +
    scale_x_continuous(labels = scales::comma) +
    scale_alpha_manual(values = c(1,0.2,1)) +
    scale_color_SU() +
    labs(x = "Population",
         y = "Adjusted admission rate per 1,000",
         title = "Variation in admission rates by ICB and underlying population",
         subtitle = paste0("Subcohort: ", subtitle_text, " | 2023"),
         colour = "",
         alpha = ""
    ) 
}

# 2x2 plots to compare admission rates in 2 cohorts on 1 graph

compare_2_2_plots <- function(cohort_1, cohort_2, cohort_1_axis, cohort_2_axis) {
  
  data_cohort_1 <-
    adjusted_rate_sub_cohorts |> 
    filter(id_clean == {{cohort_1}}, 
           year == 2023) 
  
  mean_rate_cohort_1 <- mean(data_cohort_1$adjusted_spell_rate)
  sd_rate_cohort_1 <- sd(data_cohort_1$adjusted_spell_rate)
  
  data_cohort_2 <-
    adjusted_rate_sub_cohorts |> 
    filter(id_clean == {{cohort_2}}, 
           year == 2023) 
  
  mean_rate_cohort_2 <- mean(data_cohort_2$adjusted_spell_rate)
  sd_rate_cohort_2 <- sd(data_cohort_2$adjusted_spell_rate)
  
  # Plot
  plot <-
    data_cohort_1 |> 
    mutate(std_from_mean_cohort_1 = (adjusted_spell_rate - mean_rate_cohort_1) / sd_rate_cohort_1) |> 
    select(1:2, adjusted_spell_rate, std_from_mean_cohort_1) |> 
    rename(cohort_1_adj_rate = adjusted_spell_rate) |> 
    left_join(data_cohort_2 |> 
                mutate(std_from_mean_cohort_2 = (adjusted_spell_rate - mean_rate_cohort_2) / sd_rate_cohort_2) |> 
                select(1:2, adjusted_spell_rate, std_from_mean_cohort_2) |> 
                rename(cohort_2_adj_rate = adjusted_spell_rate),
              by = c("icb23cd", "icb_name_short")) |> 
    rename(ICB = icb_name_short) |> 
    
    mutate(fill_text = 
             case_when(
               (std_from_mean_cohort_1 <= 1 & std_from_mean_cohort_1 >= -1) &
                 (std_from_mean_cohort_2 <= 1 & std_from_mean_cohort_2 >= -1) ~ "A. Within 1 SD",
               (std_from_mean_cohort_1 > 1 & std_from_mean_cohort_2 > 1) ~ "B. Above 1 SD",
               (std_from_mean_cohort_1 < -1 & std_from_mean_cohort_2 < -1) ~ "C. Below 1 SD",
               TRUE ~ "D. +/- 1 SD on one axis"
             )) |> 
    
    ggplot(aes(x = cohort_1_adj_rate, y = cohort_2_adj_rate, 
               colour = fill_text, alpha = fill_text, 
               label =  ICB
               )) +
    geom_point(size = 3) +
    
    # Cohort 1 lines
    geom_vline(xintercept = mean_rate_cohort_1, color = "blue", linetype = "dashed", linewidth = 0.3) +
    geom_vline(xintercept = mean_rate_cohort_1 + sd_rate_cohort_1, color = "red", linetype = "dotted", linewidth = 0.25) +
    geom_vline(xintercept = mean_rate_cohort_1 - sd_rate_cohort_1, color = "red", linetype = "dotted", linewidth = 0.25) +
    
    # Cohort 2 line
    geom_hline(yintercept = mean_rate_cohort_2, color = "blue", linetype = "dashed", linewidth = 0.3) +
    geom_hline(yintercept = mean_rate_cohort_2 + sd_rate_cohort_2, color = "red", linetype = "dotted", linewidth = 0.25) +
    geom_hline(yintercept = mean_rate_cohort_2 - sd_rate_cohort_2, color = "red", linetype = "dotted", linewidth = 0.25) +
    
    scale_x_continuous(labels = scales::comma) +
    scale_alpha_manual(values = c(0.3, 1, 1, 1)) + # Only shade out those within 1 SD on both axis
    scale_color_SU() +
    theme_SU() +
    labs(x = paste0("Admission rate: ", cohort_1_axis),
         y = paste0("Admission rate: ", cohort_2_axis),
         colour = "",
         alpha = ""
         )
  
  #plot
  
  ggplotly(plot, tooltip = "label") |>
    layout(
      title = list(
        text = "<b>Variation in admission rates by ICB</b><br>",
        x = 0.05,  # Center the title
        y = 0.95,
        xanchor = "left",
        font = list(size = 14)  # Title font size
      ),
      annotations = list(
        list(
          text = "<i>Comparison of adjusted admission rates by sub-cohort</i>",
          x = -0.01,  # Center the title
          y = 0.98,
          xanchor = "left",
          xref = "paper",
          yref = "paper",
          showarrow = FALSE,
          font = list(size = 12)  # Subtitle font size
          )
        )
      )
}



# Overlap in cohorts ----
overlap_data <-
  aggregate_data %>%
  filter(year(month) == 2023) |> 
  mutate(flag_falls = 
           case_when(
             !is.na(flag_falls_exp) ~ flag_falls_exp,
             !is.na(flag_fall_imp_frac) ~ flag_fall_imp_frac,
             !is.na(flag_fall_imp_tend) ~ flag_fall_imp_tend,
             TRUE ~ NA
           )) |> 
  group_by(flag_frail, flag_elderly_emergency, flag_eol, flag_falls) |> 
  summarise(spells = sum(spells)) |> 
  mutate(
    flag_frail = !is.na(flag_frail),
    flag_eol = !is.na(flag_eol),
    flag_falls = !is.na(flag_falls),
    flag_elderly_emergency = !is.na(flag_elderly_emergency)
    ) # Convert flags to logical values

# Calculate the number of spells for each combination
overlap_combinations <- 
  overlap_data %>%
  group_by(flag_frail, flag_elderly_emergency, flag_eol, flag_falls) %>%
  summarise(spells = sum(spells), .groups = 'drop')

expanded_data <- 
  overlap_combinations %>%
  rowwise() %>%
  do(data.frame(flag_frail = rep(.$flag_frail, .$spells),
                flag_eol = rep(.$flag_eol, .$spells),
                flag_falls = rep(.$flag_falls, .$spells),
                flag_elderly_emergency = rep(.$flag_elderly_emergency, .$spells)
                ))

# Convert to matrix
expanded_matrix <- as.matrix(expanded_data)

# Create a list for Venn diagram
venn_data <- list(
  Frail = which(expanded_matrix[, "flag_frail"]),
  EOL = which(expanded_matrix[, "flag_eol"]),
  Falls = which(expanded_matrix[, "flag_falls"]),
  `Emergency elderly` = which(expanded_matrix[, "flag_elderly_emergency"])
)


# ICB spell denominators ----
icb_spell_denominators <- 
  read_csv("icb_spell_denominators.csv") |>
  clean_names()

spell_proportion <- function(data, label_input) {
  
  data |> 
    group_by(der_activity_month) |> 
    summarise(cohort_spells = sum(spells)) |>  
    left_join(icb_spell_denominators |> 
                group_by(der_activity_month) |> 
                summarise(denom_spells = sum(spells)), 
              by = c("der_activity_month")) |> 
    mutate(prop = round(cohort_spells/denom_spells * 100, 3)) |>
    mutate(month = 
             as.Date(
               paste0(str_sub(der_activity_month,1,4),
                      "-",
                      str_sub(der_activity_month,5,6),
                      "-01")
             )
    ) |> 
    mutate(id = label_input) |> 
    select(month, cohort_spells, denom_spells, prop, id)
  }

# Patient characteristics ----

patchwork_function <- function(data, title_text) {
  
  # Plot 1
  plot_1 <-
    data %>% 
    drop_na(age_range) %>% 
    filter(year == 2023,
           !age_range %in% c("100-104","105-109")) %>% 
    group_by(age_range, sex) %>% 
    summarise(pts = sum(person_n)) |> 
    mutate(pts_2 = 
             case_when(sex == 1 ~ 0-pts,
                       TRUE ~ pts)) |> 
    mutate(sex = 
             case_when(sex == 1 ~ "Male",
                       TRUE ~ "Female")) %>%
    
    ggplot(aes(x = pts_2, y = age_range, fill = sex)) +
    geom_col() +
    geom_vline(xintercept = 0) +
    scale_fill_SU() +
    scale_x_continuous(labels = function(x) scales::comma(abs(x))) +
    theme(legend.position = "bottom",
          plot.subtitle = element_text(hjust = 0.5)
          ) +
    labs(x = "Patients",
         y = "Age-range",
         fill = "Sex:",
         subtitle = "Age and sex distribution",
         )
  
  
  # Plot 2
  ethnicity_imd <-
    data |> 
    drop_na(age_range) %>% 
    filter(year == 2023,
           !age_range %in% c("100-104","105-109")) %>% 
    mutate(ethnic_group = str_sub(ethnic_group, 1,1)) %>% 
    left_join(ethnicity_lookup, by = c("ethnic_group" = "Code")) |> 
    group_by(Description, imd_decile) |> 
    summarise(person_n = sum(person_n)) |> 
    ungroup()
  
  # Reshape the data for the heatmap
  heatmap_data <- dcast(ethnicity_imd, Description  ~ imd_decile, value.var = "person_n")
  
  # Convert the data to a matrix
  heatmap_matrix <- as.matrix(heatmap_data[,-1])
  rownames(heatmap_matrix) <- heatmap_data$ethnic_group
  
  plot_2 <-
    ethnicity_imd |> 
    group_by(Description) |> 
    drop_na(Description) |> 
    mutate(prop = person_n/sum(person_n) * 100) |> 
    
    ggplot(aes(x = factor(imd_decile), y = Description, fill = prop)) +
    geom_tile(alpha = 0.9) +
    scale_fill_gradient(low = "#686f73", high = "#f9bf07") +
    theme(legend.position = "bottom",
          plot.subtitle = element_text(hjust = 0.5)
    ) +
    labs(x = "IMD Decile",
         y = "Ethnic Group",
         subtitle = "IMD distribution in patient ethnicity groups",
         fill = "Proportion (%)"
    ) 
  
  plot_1 + plot_2 +
    plot_annotation(title = paste0("Sub-cohort demographics - ", title_text),
                    subtitle = "Mitigable SUS admissions | 2023")
  
  
}

# Top diagnoses by cohort
read_diagnosis_list <- function(cohort) {

  data <- 
    read_csv(paste0("top_diagnoses/", cohort, ".csv")) |> 
    clean_names() |> 
    select(4,2) |> 
    rename(Admissions = 2) |> 
    arrange(desc(Admissions)) |> 
    mutate(Proportion = round(Admissions/sum(Admissions)*100, 2)) |> 
    rename(`Primary diagnosis` = icd10_l4_desc) |> 
    head(10) |> 
    mutate(Admissions = scales::comma(Admissions))
  
  data
}

# Top diagnoses by cohort
read_procedures_list <- function(cohort) {
  
  data <- 
    read_csv(paste0("top_procedures/", cohort, ".csv")) |> 
    clean_names() |> 
    select(4,2) |> 
    rename(Admissions = 2) |> 
    arrange(desc(Admissions)) |> 
    mutate(Proportion = round(Admissions/sum(Admissions)*100, 2)) |> 
    rename(`Primary procedure` = opcs_l4_desc) |> 
    head(10) |> 
    mutate(Admissions = scales::comma(Admissions))
  
  data
}

```

# Project specifications

::: {.callout-note icon="false"}
## Ask

Size and describe patient cohorts and activity in the inpatient acute settings using National NHP model mitigators.

Generate descriptive analysis to support onward analysis by NHSE analysts covering:

-   Number of patients,
-   Levels of activity,
-   Variation across the country,
-   Characteristics of patients
:::

# Outline

::: panel-tabset
## 1. Patient characteristics

Describe the underlying patient demographics within out mitigable-activity sub-cohorts by:

-   Age and sex
-   Ethnicity and deprivation

## 2. Size and describe

Present activity totals and time-series trends for the combined mitigable-activity cohort and mitigator sub-cohorts.

Measures include:

-   Counts of patients, admissions and bed days:
    -   In combined mitigable-activity cohort
    -   By mitigator sub-cohorts
    -   By ICB
-   Cohort activity as a share of all admissions

## 3. Use variation

Explore variation in age- and sex-adjusted admission rates to:

-   Examine the range of activity by ICB
-   Assess whether systems in place are particularly good or bad at treating mitigable activity

**Additional analysis:** Incorporate underlying population need and/or disease prevalence (weighting by deprivation or regression controlling for population variables)

Apply survival analysis techniques to assess differences in care outcomes by patient groups and ICBs. Care outcomes include:

-   Readmission post-acute inpatient care
-   Motality post-acute inpatient care

**Additional analysis:**

-   Summary of care services/types that patients from mitigable cohorts are in contact with in the year prior to death.
-   Comparison of location of death in mitigable patients that recieved acute inpatient care in the year prior to death.

## 4. Mitigator overlap

Display the overlap of patients in mitigable cohorts by combined and sub-cohorts.

Outputs include:

-   Underlying SQL queries
-   Code
-   Data extract
-   Quarto report output

## 5. Additional deliverable(s)

Ranked table of activity in each mitigator with activity counts and proportion of activity
:::

<hr style="border-top: 2px solid #f9bf07">

# Cohort definition

In keeping with the wider context of this work, we have recreated mitigable activity cohorts used elsewhere in the New Hospital Programme (NHP) Demand and Capacity Model. Our extracts differ primarily in datasource, where the NHP model is built using Hosptial Episode Data (HES) while we have accessed Secondary Uses Service (SUS) data via the National Commissioning Data Repository (NCDR) portal. As such, we have included data between April 2018 and August 2024 and include activity from across England.

Full SQL queries used to define our patient cohorts can be found on the project github repository: [sql_queries](https://github.com/The-Strategy-Unit/Community_Strategies/tree/main/sql).

Details of the original NHP queries can be found here: [NHP mitigators](https://connect.strategyunitwm.nhs.uk/nhp/project_information/modelling_methodology/activity_mitigators/inpatient_activity_mitigators.html).

::: panel-tabset

### Emergency elderly

Admissions included in our emergency elderly cohort are all emergency admissions in patients aged 75 and older (Admission_Method LIKE "2%" - details [here](https://www.datadictionary.nhs.uk/attributes/admission_method.html)). 


### Frail

Admissions included in our Frail cohort are all emergency admissions in patients aged 65 and older that include a Spell_Primary_Diagnosis code indicating frailty as the cause of the admission. This is different to the NHP frailty definition which focuses only on patients aged 75 years and older.

::: {.callout-caution collapse="true"}
## Expand To view ICD10 frailty diagnosis codes:

LEFT([Spell_Primary_Diagnosis],3) IN (
'F00','G81','G30','I69','R29','R39','F05','W19',
'S00','R31','B96','R41','R26','I67','R56','R40',
'T83','S06','S42','E87','M25','E86','R54','Z50',
'F03','W18','Z75','F01','S80','L03','H54','E53',
'Z60','G20','R55','S22','K59','N17','L89','Z22',	
'B95','L97','R44','K26','I95','N19','A41','Z87',	
'J96','X59','M19','G40','M81','S72','S32','E16',	
'R94','N18','R33','R69','N28','R32','G31','Y95',	
'S09','R45','G45','Z74','M79','W06','S01','A04',	
'A09','J18','J69','R47','E55','Z93','R02','R63',	
'H91','W10','W01','E05','M41','R13','Z99','U80',	
'M80','K92','I63','N20','F10','Y84','R00','J22',	
'Z73','R79','Z91','S51','F32','M48','E83','M15',	
'D64','L08','R11','K52','R50'
) 

:::

### Falls 

Our falls cohort included any patient aged 65 years or older that had an emergency admission and a primary diagnosis indicating a fall (expliticly or implicitly).

::: {.callout-caution collapse="true"}
## Expand To view ICD10 falls diagnosis codes:

**Explicit**: 
Spell_Primary_Diagnosis LIKE 'W[01][0-9]%'

**Implicit - fractures**: 
     (
          Spell_Primary_Diagnosis LIKE 'M48[45]%'
        OR
          Spell_Primary_Diagnosis LIKE 'M80[0-589]%'
        OR
          Spell_Primary_Diagnosis LIKE 'S22[01]%'
        OR
          Spell_Primary_Diagnosis LIKE 'S32[0-47]%'
        OR
          Spell_Primary_Diagnosis LIKE 'S42[234]%'
        OR
          Spell_Primary_Diagnosis LIKE 'S52%'
        OR
          Spell_Primary_Diagnosis LIKE 'S620%'
        OR
          Spell_Primary_Diagnosis LIKE 'S72[0-48]%'
        OR
          Spell_Primary_Diagnosis LIKE 'T08X%'
      ) 
	  AND
	  Spell_Primary_Diagnosis NOT LIKE '[VWXY]%'
	  
**Implicit - Tendency to fall**: 
Spell_Primary_Diagnosis LIKE 'R296%'
:::


### End of life

End of life spells were identified by filtering admissions on Dischare_Method = 4 (Outcome - Patient died. Details  [here](https://www.datadictionary.nhs.uk/data_elements/discharge_method_code__hospital_provider_spell_.html)) where the length of stay (Der_Spell_LoS) was less than 14 days and no procedure was undertaken (Der_Procedure_All IS NULL). 



:::

<hr style="border-top: 2px solid #f9bf07">

# Patient characteristics

## Demographics

In both our emergency elderly and frail patient groups, it is apparent that within the Black and Asian ethnic groups, the majority of patients are clustered in the most deprived IMD deciles (particularly deciles 2-4). However, when considering emergency elderly admissions, in larger group of White British patients, the proportion of patients gradually increases with affluence.

::: panel-tabset
### Emergency elderly

```{r}
patchwork_function(aggregate_data_elderly_emergency, "Emergency elderly patients")
```

### Frail

```{r}
patchwork_function(aggregate_data_frail, "Frail patients")
```

### Falls

```{r}
patchwork_function(aggregate_data_falls, "Falls patients")
```

### End-of-life

```{r}
patchwork_function(aggregate_data_eol, "End-of-life patients")
```
:::

## Most frequent diagnoses

::: panel-tabset
### Emergency elderly

```{r}
create_dt(read_diagnosis_list("frail"))
```

### Frail

```{r}
create_dt(read_diagnosis_list("emergency_elderly"))
```

### Falls

```{r}
create_dt(read_diagnosis_list("falls"))
```

### End-of-life

```{r}
create_dt(read_diagnosis_list("eol"))
```
:::

## Most frequent procedures

::: panel-tabset
### Emergency elderly

```{r}
create_dt(read_procedures_list("frail"))
```

### Frail

```{r}
create_dt(read_procedures_list("emergency_elderly"))
```

### Falls

```{r}
create_dt(read_procedures_list("falls"))
```

### End-of-life

```{r}
create_dt(read_procedures_list("eol"))
```
:::

<hr style="border-top: 2px solid #f9bf07">

# Size and describe

The following section will describe mitigable-activity in our combined cohort and mitigator-specific sub-cohorts.

We have created our data set by applying existing strategies developed under the New Hospitals Program (NHP) that filter patient-level hospital activity data (Secondary Uses Service data - SUS) which center on patient groups and pathways which might reasonably be suitable for treatment in the community.

The following patient groups have been identified with such mitigation strategies:

-   Admissions in frail elderly patients (65 years +)
-   Emergency admissions in over 75's
-   End-of-life admissions: spells that ended with the patient dying and lasted less that 14 days
-   Admissions post-falls (slips, trips and falls)

## Combined mitigable activity cohort

::: panel-tabset
### Time series

```{r}
aggregate_data_core_cohorts %>%  
  group_by(month) %>% 
  summarise(`1. Spells` = sum(spells),
            `2. Individuals` = sum(person_n),
            `3. Bed days` = sum(los_sum)) %>% 
  pivot_longer(cols = -month) %>% 
  
  ggplot(aes(x = month, y = value, colour = name)) +
  #geom_point() +
  geom_line(linewidth = 1) +
  facet_wrap(~name, scales = "free_y", nrow = (3)) +
  scale_y_continuous(labels = scales::comma) +
  scale_color_SU() +
  theme(legend.position = "none",
        strip.background = element_rect(fill = NA, colour = "grey"),
        axis.title.y = element_blank()
        ) +
  labs(x = "Month",
       title = "National trend in mitigable acute inpatient admissions",
       subtitle = "SUS Apr 2018 - Aug 2024")
```

### Annual table

```{r}
create_dt(
aggregate_data_core_cohorts %>%  
  mutate(year  = lubridate::year(month)) %>% 
  group_by(year) %>% 
  summarise(`1. Spells` = sum(spells),
            `2. Individuals` = sum(person_n),
            `3. Bed days` = sum(los_sum)
            ) %>%
  pivot_longer(cols = -year) %>% 
  pivot_wider(id_cols = year, names_from = name, values_from = value)
)
```
:::

Activity in the combined data set has gradually returned to pre-pandemic levels. While the increase in post-pandemic growth in spells has been gradual between 2021 and 2024, the length of time a person from our cohort has stayed in hospital has risen quickly to surpass pre-pandemic levels. The graph above suggests a change in the average length of stay particularly between 2021 and 2023, after which bed-days associated with our activity have plateaued.

Is the change in bed days linked to an easing of pandemic-related practices around length of treatment and/or discharge?

::: panel-tabset
### Share of total admissions

```{r}
aggregate_data_core_cohorts |> 
  group_by(der_activity_month) |> 
  summarise(cohort_spells = sum(spells)) |>  
  left_join(icb_spell_denominators |> 
              group_by(der_activity_month) |> 
              summarise(denom_spells = sum(spells)), 
            by = c("der_activity_month")) |> 
  mutate(prop = round(cohort_spells/denom_spells * 100, 3)) |>
  mutate(month = 
           as.Date(
             paste0(str_sub(der_activity_month,1,4),
                    "-",
                    str_sub(der_activity_month,5,6),
                    "-01")
             )
         ) |> 
  
  ggplot(aes(x = month, y = prop)) +
  #geom_point() +
  geom_line() +
  geom_smooth(method = "loess", colour = "#f9bf07") +
  labs(y = "Cohort proportion",
       title = "Cohort spells as a proportion of total spells",
       subtitle = "NHS England admissions denominator | 2018-24")
```

### Data

```{r}
create_dt(
  aggregate_data_core_cohorts |> 
    group_by(der_activity_month) |> 
    summarise(cohort_spells = sum(spells)) |>  
    left_join(icb_spell_denominators |> 
                group_by(der_activity_month) |> 
                summarise(denom_spells = sum(spells)), 
              by = c("der_activity_month")) |> 
    mutate(prop = cohort_spells/denom_spells * 100) |>
    mutate(month = 
             as.Date(
               paste0(str_sub(der_activity_month,1,4),
                      "-",
                      str_sub(der_activity_month,5,6),
                      "-01")
               )
           ) |> 
    select(month, cohort_spells, denom_spells, prop)
)
```
:::

We contextualise the activity identified in our cohort alongside wider NHS-funded healthcare delivery. The combined sum of admissions from our emergency eldery, frail, falls and end-of-life cohorts accounts to approx 11% of total admissions. The trend varied most significantly during the COVID-19 pandemic and has since returned to pre-pandemic levels, allowing for consistent seasonal peaked during the winter months.

## Sub-cohorts

By splitting our activity according to mitigation strategy, we can compare trends in patient group. The volume of activity by cohort differs significantly with Emergency Elderly and Frail activity accounting for the larger shares of total admissions and bed days.

::: panel-tabset
### Fixed axis

```{r}
aggregate_data_frail %>% 
  group_by(month) %>% 
  summarise(frail_spells = sum(spells),
            frail_ind = sum(person_n),
            frail_bed_days = sum(los_sum)) %>% 
  left_join(
    aggregate_data_eol %>% 
      group_by(month) %>% 
      summarise(eol_short_spells = sum(spells),
                eol_short_ind = sum(person_n),
                eol_short_bed_days = sum(los_sum)),
    by = "month"
    ) %>%
  left_join(
    aggregate_data_falls %>% 
      group_by(month) %>% 
      summarise(falls_spells = sum(spells),
                falls_ind = sum(person_n),
                falls_bed_days = sum(los_sum)),
    by = "month"
  ) %>%
  left_join(
    aggregate_data_elderly_emergency %>% 
      group_by(month) %>% 
      summarise(emergency_elderly_spells = sum(spells),
                emergency_elderly_ind = sum(person_n),
                emergency_elderly_bed_days = sum(los_sum)),
    by = "month"
  ) %>%
  #left_join(
  #  aggregate_data_amb_acute %>% 
  #    group_by(month) %>% 
  #    summarise(amb_acute_spells = sum(spells),
  #              amb_acute_bed_days = sum(los_sum)),
  #  by = "month"
  #) %>%
  #left_join(
  #  aggregate_data_amb_chronic %>% 
  #    group_by(month) %>% 
  #    summarise(amb_chronic_spells = sum(spells),
  #              amb_chronic_bed_days = sum(los_sum)),
  #  by = "month"
  #) %>%
  #left_join(
  #  aggregate_data_amb_vacc_prev %>% 
  #    group_by(month) %>% 
  #    summarise(amb_vaccine_spells = sum(spells),
  #              amb_vaccine_bed_days = sum(los_sum)),
  #  by = "month"
  #) %>%
  #left_join(
  #  aggregate_data_eol_1_year %>% 
  #    group_by(month) %>% 
  #    summarise(eol_1_year_spells = sum(spells),
  #              eol_1_year_bed_days = sum(los_sum)),
  #  by = "month"
  #) %>%
  
  pivot_longer(cols = -month) %>% 
  filter(month != as.Date("2024-10-01")) |> 
  mutate(type = case_when(str_detect(name, "bed_days") ~ "3. Bed days",
                          str_detect(name, "ind") ~ "2. Individuals",
                          TRUE ~ "1. Spells"),
         Cohort = case_when(str_detect(name, "frail") ~ "1. Frail",
                            str_detect(name, "eol_short") ~ "3. End of life",
                            str_detect(name, "falls") ~ "4. Falls",
                            str_detect(name, "emergency_elderly") ~ "2. Emergency elderly",
                            str_detect(name, "amb_acute") ~ "5. Ambulatory - Acute",
                            str_detect(name, "amb_chronic") ~ "6. Ambulatory - Chronic",
                            str_detect(name, "amb_vacc") ~ "7. Ambulatory - Vaccine preventable",
                            str_detect(name, "eol_1_year") ~ "8. End of life - 1 year"
                            )
         ) %>% 
  
  ggplot(aes(x = month, y = value, colour = Cohort)) +
  #geom_point() +
  geom_line(linewidth = 1) +
  facet_wrap(~type, scales = "free_y") +
  scale_y_continuous(labels = scales::comma) +
  scale_color_SU() +
  theme(axis.text.x = element_text(angle = 90),
        strip.background = element_rect(fill = NA, colour = "grey"),
        axis.title.y = element_blank()
        ) +
  labs(x = "Month",
       title = "National trend in acute inpatient admissions",
       subtitle = "SUS Apr 2018 - Aug 2024")
```

### Free axis

```{r}
aggregate_data_frail %>% 
  group_by(month) %>% 
  summarise(frail_spells = sum(spells),
            frail_ind = sum(person_n),
            frail_bed_days = sum(los_sum)) %>% 
  left_join(
    aggregate_data_eol %>% 
      group_by(month) %>% 
      summarise(eol_short_spells = sum(spells),
                eol_short_ind = sum(person_n),
                eol_short_bed_days = sum(los_sum)),
    by = "month"
  ) %>%
  left_join(
    aggregate_data_falls %>% 
      group_by(month) %>% 
      summarise(falls_spells = sum(spells),
                falls_ind = sum(person_n),
                falls_bed_days = sum(los_sum)),
    by = "month"
  ) %>%
  left_join(
    aggregate_data_elderly_emergency %>% 
      group_by(month) %>% 
      summarise(emergency_elderly_spells = sum(spells),
                emergency_elderly_ind = sum(person_n),
                emergency_elderly_bed_days = sum(los_sum)),
    by = "month"
  ) %>%
  pivot_longer(cols = -month) %>% 
  
  mutate(type = case_when(str_detect(name, paste(c("spells", "ind"), collapse = "|")) ~ "1. Spells & Individuals",
                          #str_detect(name, "bed_days") ~ "2. Bed days",
                          TRUE ~ "2. Bed days"),
         Cohort = case_when(str_detect(name, "frail") ~ "Frail",
                            str_detect(name, "eol") ~ "End of life",
                            str_detect(name, "falls") ~ "Falls",
                            str_detect(name, "emergency") ~ "Emergency elderly"),
         label = case_when(str_detect(name, "spells") ~ "1. Spells",
                           str_detect(name, "ind") ~ "2. Individuals",
                           str_detect(name, "bed_days") ~ "3. Bed days")
         ) %>% 
  #mutate(Cohort = factor(Cohort, levels = "Frail", "Emergency elderly", "End of life", "Falls")) %>% 
  
  ggplot(aes(x = month, y = value, colour = label, group = name)) +
  #geom_point() +
  geom_line(linewidth = 1) +
  facet_wrap(Cohort~type, scales = "free", ncol = 2) +
  scale_y_continuous(labels = scales::comma) +
  scale_color_SU() +
  theme(axis.text.x = element_text(angle = 90),
        strip.background = element_rect(fill = NA, colour = "grey"),
        axis.title.y = element_blank()
        #legend.position = "none"
        ) +
  labs(x = "Month",
       colour = "",
       title = "National trend in acute inpatient admissions",
       subtitle = "SUS Apr 2018 - Aug 2024")
```

### Table

```{r}
create_dt(
  aggregate_data_frail %>% 
  group_by(year) %>% 
  summarise(frail_spells = sum(spells),
            frail_ind = sum(person_n),
            frail_bed_days = sum(los_sum)) %>% 
  left_join(
    aggregate_data_eol %>% 
      group_by(year) %>% 
      summarise(eol_spells = sum(spells),
                eol_ind = sum(person_n),
                eol_bed_days = sum(los_sum)),
    by = "year"
  ) %>%
  left_join(
    aggregate_data_falls %>% 
      group_by(year) %>% 
      summarise(falls_spells = sum(spells),
                falls_ind = sum(person_n),
                falls_bed_days = sum(los_sum)),
    by = "year"
  ) %>%
  left_join(
    aggregate_data_elderly_emergency %>% 
      group_by(year) %>% 
      summarise(emergency_elderly_spells = sum(spells),
                emergency_elderly_ind = sum(person_n),
                emergency_elderly_bed_days = sum(los_sum)),
    by = "year"
  ) %>%
  pivot_longer(cols = -year) %>%
  mutate(value = scales::comma(value)) |> 
  pivot_wider(id_cols = year, names_from = name, values_from = value) |> 
  rename(
    `Frail - spells` = frail_spells,
    `Frail - individuals` = frail_ind,
    `Frail - bed days` = frail_bed_days,
    `EoL - spells` = eol_spells,
    `EoL - individuals` = eol_ind,
    `Eol - bed days` = eol_bed_days, 
    `Falls - spells` = falls_spells,
    `Falls - individuals` = falls_ind,
    `Falls - bed days` = falls_bed_days,
    `Emergency elderly - spells` = emergency_elderly_spells,
    `Emergency elderly - individuals` = emergency_elderly_ind,
    `Emergency elder - bed days` = emergency_elderly_bed_days
  )
)
```
:::

Though the post-pandemic trend in frail activity is growing at a steady rate (from c.90,000 per month in 2022 up to c.100,000 per month since), there has been a recent stepped increase in emergency elderly admissions in the last year (stable around 140,000 per month between 2021-23 but increased to 160,000 since start of 2024). For both of these cohorts, the length of stay associated with these admissions has grown substantially and is higher than pre-pandemic levels in the emergency elderly cohort.

Similarly, the length of stay associated with falls admissions has surged while the underlying activity has remained stable.

Excluding pandemic related surges in end of life care in secondary settings, trends in activity and bed days are stable and correlated.

::: panel-tabset
### Share of total admissions

```{r}
spell_proportion(aggregate_data_elderly_emergency, "1. Emergency elderly") |> 
  union_all(spell_proportion(aggregate_data_frail, "2. Frail")) |> 
  union_all(spell_proportion(aggregate_data_falls, "3. Falls")) |>
  union_all(spell_proportion(aggregate_data_eol, "4. End of life")) |> 
  
  #union_all(spell_proportion(aggregate_data_amb_acute, "5. Ambulatory - acute")) |> 
  #union_all(spell_proportion(aggregate_data_amb_chronic, "6. Ambulatory - chronic")) |> 
  #union_all(spell_proportion(aggregate_data_amb_vacc_prev, "7. Ambulatory - vaccine preventable")) |> 
  #union_all(spell_proportion(aggregate_data_eol_1_year, "8. End of life - 1 year")) |> 
  
  ggplot(aes(x = month, y = prop, colour = id)) +
  geom_point() +
  geom_line() +
  facet_wrap(~id#, scales = "free_y"
  ) +
  geom_smooth(method = "loess") +
  scale_color_SU() +
  theme(strip.background = element_rect(fill = NA, colour = "grey"),
        axis.title.x = element_blank(),
        legend.position = "none"
  ) +
  labs(y = "Cohort proportion",
       colour = "Cohort:",
       title = "Spells as a proportion of total spells",
       subtitle = "NHS England | 2018-24")
```

### Free scale

```{r}
spell_proportion(aggregate_data_elderly_emergency, "1. Emergency elderly") |> 
  union_all(spell_proportion(aggregate_data_frail, "2. Frail")) |> 
  union_all(spell_proportion(aggregate_data_falls, "3. Falls")) |>
  union_all(spell_proportion(aggregate_data_eol, "4. End of life")) |> 
  
  #union_all(spell_proportion(aggregate_data_amb_acute, "5. Ambulatory - acute")) |> 
  #union_all(spell_proportion(aggregate_data_amb_chronic, "6. Ambulatory - chronic")) |> 
  #union_all(spell_proportion(aggregate_data_amb_vacc_prev, "7. Ambulatory - vaccine preventable")) |> 
  #union_all(spell_proportion(aggregate_data_eol_1_year, "8. End of life - 1 year")) |> 
  
  ggplot(aes(x = month, y = prop, colour = id)) +
  geom_point() +
  geom_line() +
  facet_wrap(~id, scales = "free_y"
  ) +
  geom_smooth(method = "loess") +
  scale_color_SU() +
  theme(strip.background = element_rect(fill = NA, colour = "grey"),
        axis.title.x = element_blank(),
        legend.position = "none"
  ) +
  labs(y = "Cohort proportion",
       colour = "Cohort:",
       title = "Spells as a proportion of total spells",
       subtitle = "NHS England | 2018-24")
```

### Table

```{r}
create_dt(
  spell_proportion(aggregate_data_elderly_emergency, "1. Emergency elderly") |> 
    union_all(spell_proportion(aggregate_data_frail, "2. Frail")) |> 
    union_all(spell_proportion(aggregate_data_falls, "3. Falls")) |>
    union_all(spell_proportion(aggregate_data_eol, "4. End of life"))
)
```
:::

When we separate our patient cohorts, we're reminded of the comparative differences in scale. The emergency elderly and frail cohorts account for approximately 10 and 6.5 per cent of total inpatient hospital activity respectively, while falls and end of life care are both below 1% total admissions. All cohorts saw increases in proportion during the pandemic as activity in other patient groups reduced.

## ICB trends

```{r}

icb_23_shp %>% 
  select(icb23cd) %>% 
  left_join(
    aggregate_data_core_cohorts %>%  
      mutate(year = year(month)) %>% 
      filter(year == 2023) %>% 
      group_by(icb_name_short) %>% 
      summarise(spells = sum(spells)) %>% 
      mutate(icb_code = str_extract(icb_name_short, "^[^:]+")) %>% 
      left_join(icb_lookup, by = c("icb_code" = "icb23cdh")),
    by = "icb23cd"
    ) %>% 

# Plot chloropleth map 
ggplot() +
  geom_sf(aes(fill = spells)) +
  scale_fill_gradient(low = "yellow", high = "red", 
                      name = "Admissions",
                      labels = function(x) format(x, big.mark = ",", scientific = FALSE)) +
  theme(panel.grid = element_blank(),
        axis.text = element_blank()
        ) +
  labs(title = "Mitigable activity by ICB",
       subtitle = "SUS 2023 | Admission identified: Frail, emergency elderly, falls & end of life")

```

In absolute volume terms, ICB's differ significantly in the number of mitigable activity undertaken in secondary care settings though this may largely be a function of differences in underlying population distribution of elderly people.

::: panel-tabset
### Volumes - facet

```{r}
aggregate_data_core_cohorts %>%  
  group_by(month, icb_name_short) %>% 
  summarise(spells = sum(spells)) %>% 
  filter(icb_name_short != "NULL") %>% 
  mutate(icb_name_short = str_extract(icb_name_short, "(?<=NHS).*")) %>% 
  mutate(icb_name_short = str_remove_all(icb_name_short, " ICB")) |> 
  
  ggplot(aes(x = month, y = spells, group = icb_name_short)) +
  #geom_point() +
  geom_line() +
  facet_wrap(~str_wrap(icb_name_short, 20), scales = "free_y") +
  scale_y_continuous(labels = scales::comma) +
  labs(x = "Month",
       y = "Spells",
       title = "ICB trends in mitigable acute inpatient admissions",
       subtitle = "SUS Apr 2018 - Aug 2024")
```

### Line

```{r}
aggregate_data_core_cohorts %>%  
  group_by(month, icb_name_short) %>% 
  summarise(spells = sum(spells)) %>% 
  filter(icb_name_short != "NULL") %>% 
  mutate(icb_name_short = str_extract(icb_name_short, "(?<=NHS).*")) %>% 
  
  ggplot(aes(x = month, y = spells)) +
  #geom_point() +
  geom_line(aes(group = icb_name_short), colour = "#2c2825") +
  geom_smooth(aes(group = 1), method = "loess", colour = "#f9bf07") +
  scale_y_continuous(labels = scales::comma) +
  labs(x = "Month",
       y = "Spells",
       title = "ICB trends in mitigable acute inpatient admissions",
       subtitle = "SUS Apr 2018 - Aug 2024")
```

### Table

```{r}
create_dt(
  aggregate_data_core_cohorts %>%  
  group_by(month, icb_name_short) %>% 
  summarise(spells = sum(spells)) %>% 
  filter(icb_name_short != "NULL") %>% 
  mutate(icb_name_short = str_extract(icb_name_short, "(?<=NHS).*"))
)
```
:::

The trends in mitigable activity by ICB vary significantly.

Data quality issues are apparent in a handful of ICB's:

-   Frimley,
-   Shropshire,
-   Telford and Wrekin,
-   Cambs & Peterborough and
-   Dorset.

While some ICB's display sustained reductions in mitigable admissions:

-   South West London,
-   West Yorkshire and
-   Lincolnshire.

Others show significant growth in activity: 

-   Devon,
-   Staffordshire & Stoke-on-Trent and
-   South Yorkshire for example


Are boundary changes impactful here?

### Share of total activity in ICB

Below we display the range of variation in the mitigable activity as a proportion of total admissions by ICB. Each boxplot represents monthly proportions and particularly high or low outlier ICB's are identified.

```{r}
icb_spells_proportion_total <-
  aggregate_data_core_cohorts |> 
  group_by(der_activity_month, icb_name_short) |> 
  summarise(cohort_spells = sum(spells)) |> 
  left_join(icb_spell_denominators, by = c("der_activity_month", "icb_name_short")) |> 
  mutate(prop = cohort_spells/spells * 100) |> 
  mutate(month = 
           as.Date(
             paste0(str_sub(der_activity_month,1,4),
                    "-",
                    str_sub(der_activity_month,5,6),
                    "-01")
             )
         ) |>
  mutate(icb_name_clean = str_sub(icb_name_short, 10,100)) |> 
  mutate(icb_name_clean = str_remove_all(icb_name_clean, " ICB"))

icb_spells_proportion_total |> 
  drop_na(icb_name_clean) |> 
  ggplot(aes(x = month, y = prop, group = month,)) +
  geom_boxplot(width = 5 ) +
 
  #stat_boxplot(geom = "errorbar", width = 0.2) +
  geom_line(data = icb_spells_proportion_total |> 
               filter(icb_name_clean %in% c(#"Frimley",
                                            "Dorset",
                                            "North Central London",
                                            "Sussex",
                                            "Somerset",
                                            "North East London")),
             aes(colour = icb_name_clean, group = icb_name_clean), linewidth = 1) +
  scale_color_SU() +
  theme(axis.title.x = element_blank()) +
  labs(y = "Cohort proportion",
       colour = "Outlier ICB:",
       title = "Cohort admissions as a proportion of total activity",
       subtitle = "ICB 2018-24")
```

<hr style="border-top: 2px solid #f9bf07">

# Use variation

To account for differences in the underlying number of elderly patients in each ICB, we standardised admission rates by age and sex to assess activity by ICB, if the age and sex structure of each ICB mirrored that of England (mid 2023).

::: {.callout-note icon="false"}
## Box Plot Structure

-   **Box:** The box itself represents the middle 50% of the data, with the bottom line marking the 25th percentile (Q1) and the top line marking the 75th percentile (Q3).
-   **Line inside the box:** This is the median, which divides the box into two equal halves.
-   **Whiskers:** The whiskers extend from the box to show the range of the data, excluding outliers.
-   **Dots:** Individual data points that fall outside the whiskers are considered outliers.
:::

Admission rates range from 25-55 admissions per 1,000 population between our ICB's during our data collection period. While the general trend in admission rates is increasing, there is significant variation between ICB's over the last 5 years, having accounted for age and sex structure.

::: panel-tabset
## Boxplot

```{r}
standardisation_function(aggregate_data) %>% 
  filter(year != 2024) %>% 
  ggplot(aes(x = year, y = adjusted_spell_rate, group = year)) +
  geom_point() +
  geom_boxplot() +
  labs(x = "Year",
       y = "Adjusted admission rate",
       title = "Adjusted admission rate per 1,000 population by ICB",
       subtitle = "SUS mitigable activity 2018-23 | England standard population")
```

## Line chart

```{r}
standardisation_function(aggregate_data) %>% 
  filter(year != 2024) %>% 
  ggplot(aes(x = year, y = adjusted_spell_rate)) +
  #geom_line(aes(group = icb_name_short), alpha = 0.5) +
  
  geom_smooth(aes(group = icb_name_short), span = 0.4, size = 0.5, colour = "grey") +
  
  geom_smooth(aes(group = 1), colour = "#5881c1", size = 1.2) +
  labs(x = "Year",
       y = "Adjusted admission rate",
       title = "Adjusted admission rate per 1,000 population by ICB",
       subtitle = "SUS mitigable activity 2018-23 | England standard population")
```
:::

## Admission rate by sub-cohort

The overall trends in Emergency elderly and Frail patients are similar as the frail cohort is largely a sub-set of the emergency elderly group; increasing admission rates are seen in the post-pandemic period in both. While the trend in falls admission rates is more stable, more statistical outliers exist. The trend in end of life admission rates is reducing post-pandemic but has not returned to pre-pandemic levels and considerable variation by ICB is present where it wasn't before the pandemic.

```{r}
adjusted_rate_sub_cohorts %>% 
  filter(year != 2024,
         id != "rate") %>% 
  ggplot(aes(x = year, y = adjusted_spell_rate, group = year, fill = id_clean)) +
  #geom_point() +
  geom_boxplot() +
  facet_wrap(~id_clean, scales = "free_y") +
  theme_minimal() +
  theme(legend.position = "none",
        strip.background = element_rect(fill = NA, colour = "grey")) +
  labs(y = "Adjusted admission rate",
       title = "Mitigable activity - Admission rate per 1,000 population by sub-cohort",
       subtitle = "SUS 2018-24")

```

To assess the geographic distribution of admission rates we split our data into quintiles. Quintile 1 includes the ICB's with the lowest 20% admission rates, while quintile 5 includes ICB's with the highest 20%.

Visually we can identify that ICB's in the higher admission rate quintiles for emergency elderly and frail cohorts are similar and cluster in the northern and midlands geographies. Where are ICB's with the highest falls admission rates are located around London and ICB's most often utilising secondary settings for end of life care are seen in the east, around Bristol, Worcester and Birmingham.

```{r}
combined_map_icb_rate_cohorts
```

Alternatively we can visualise admission rates in our sub-cohorts by ICB using the below bar chart.

```{r}
adjusted_rate_sub_cohorts |> 
  filter(year == 2023) |> 
  select(-icb23cd, -year, -id) |> 
  filter(id_clean != "E. Combined cohort") |> 
  group_by(id_clean) |> 
  mutate(quintile = 
           case_when(adjusted_spell_rate >= quantile(adjusted_spell_rate, 0.8)~ "5",
                     adjusted_spell_rate <= quantile(adjusted_spell_rate, 0.2)~ "1",
                     TRUE ~ "2-4")) |>
  mutate(facet = "") |>
  mutate(icb_name_clean = str_sub(icb_name_short, 10,100)) |> 
  mutate(icb_name_clean = str_remove_all(icb_name_clean, " ICB"))|> 
  mutate(icb_name_clean = fct_reorder(icb_name_clean, adjusted_spell_rate)) |>
  
  ggplot(aes(y = icb_name_clean, x = adjusted_spell_rate, fill = quintile)) +
  geom_col() +
  facet_grid(facet~id_clean, scales = "free_x") +
  scale_fill_SU() +
  theme(strip.background.x = element_rect(fill = NA, colour = "grey")) +
  labs(x = "Adjusted admission rate",
       y = "ICB",
       fill = "Quintile:",
       title = "Admission rate by ICB and sub-cohort",
       subtitle = "Age and sex adjusted rate per 1,000 population | SUS admissions 2023")
```

Or we can compare cohort activity by ICB using an average Length of Stay measure:

```{r}
avg_los_icb_cohort %>% 
  filter(icb_name_short != "NULL") %>% 
  mutate(emergency_elderly_avg_los = as.numeric(emergency_elderly_avg_los),
         falls_exp_avg_los = as.numeric(falls_exp_avg_los),
         falls_imp_frac_avg_los = as.numeric(falls_imp_frac_avg_los),
         falls_imp_tend_avg_los = as.numeric(falls_imp_tend_avg_los),
         eol_avg_los = as.numeric(eol_avg_los)
         ) %>% 
  mutate(icb_name_clean = str_sub(icb_name_short, 10,100)) |> 
  mutate(icb_name_clean = str_remove_all(icb_name_clean, " ICB"))|>
  mutate(icb_name_clean = fct_reorder(icb_name_clean, emergency_elderly_avg_los)) |>
  select(-icb_name_short) %>% 
  pivot_longer(-icb_name_clean) %>% 
  group_by(name) %>% 
  mutate(quintile = 
           case_when(value >= quantile(value, 0.8, na.rm = TRUE)~ "5",
                     value < quantile(value, 0.2, na.rm = TRUE)~ "1",
                     TRUE ~ "2-4")) %>% 
  mutate(facet = "") |>
  mutate(name = 
           case_when(
             name == "emergency_elderly_avg_los" ~ "1. Emergency elderly",
             name == "frail_avg_los" ~ "2. Frail",
             name == "eol_avg_los" ~ "3. End of life",
             name == "falls_imp_frac_avg_los" ~ "4. Falls - Implicit fracture",
             name == "falls_imp_tend_avg_los" ~ "5. Falls - Implicit tendency",
             name == "falls_exp_avg_los" ~ "6. Falls - Explicit"
           )) %>% 
  filter(name != "6. Falls - Explicit") %>% 
  
  ggplot(aes(y = icb_name_clean, x = value, fill = quintile)) +
  geom_col() +
  facet_grid(facet~str_wrap(name, 20), scales = "free_x") +
  scale_fill_SU() +
  scale_x_continuous(breaks = scales::pretty_breaks()) + # Ensure whole numbers on x-axis
  theme(strip.background.x = element_rect(fill = NA, colour = "grey")) +
  labs(x = "Average Length of Stay (days)",
       y = "ICB",
       fill = "Quintile:",
       title = "Average Length of Stay by ICB and sub-cohort",
       subtitle = "SUS admissions 2023")
```

It is important to note however, there are limitations is using length of stay as a measure by which to compare ICB's. Specifically : 

-   The range in length of stay values is small across most cohorts so the distinction between quintiles is often +/- a single day in average length of stay. 
-   The underlying data masks smaller differences in treatment length because stays are rounded to whole numbers. 
-   Similarly, presenting average length of stay does not illustrate the range of values within ICB's that may be focused on specific populations or care types. 
-   Finally, presenting an annual average figure must be considered alongside data quality issues that affect only part of 2023 (the data collection period in question for the above graph) - Frimley ICB for example, is shown to have significant data quality issues affecting the later half of 2023, which need to be taken into account when interpreting the above visualisation.

### Funnel plots

If we revert to considering admission rates per 1,000 population by patient sub-cohort and ICB, below we have identified the ICB's that sit outside (+/-) 1 standard deviation from the mean admission rate and have plotted them according to underlying population size.

For frail, emergency elderly and falls cohorts, we note that those ICB's that are demonstrating particularly low admission rates, they are often clustered in the lower left quadrant of the graph, indicating smaller populations (around 1,000,000 residents). However, the inverse is true for patients who died in hospital within 14 days of admission, where by the ICB's with higher admission rates are more likely to have smaller resident populations.

::: panel-tabset
#### Frail

```{r}
funnel_plot_function("A. Frail", "Frail")
```

#### Emergency elderly

```{r}
funnel_plot_function("B. Emergency elderly", "Emergency elderly")
```

#### Falls

```{r}
funnel_plot_function("C. Falls", "Falls")
```

#### End of life

```{r}
funnel_plot_function("D. End of life", "End of life")
```
:::

### 2x2 plots

Finally when considering variation in use of acute care, we assess ICB admission rates in our patient cohorts against other cohorts to attempt to demonstrate the trade-offs at play for the system. Click through tabs for plot description:









::: panel-tabset
#### Frail x Emergency elderly

-   **Frail x Emergency elderly**: The cohorts are highly correlated (as the frail cohort is largely a subset of the emergency elderly cohort). ICB's tend to be positioned similarly on both measures - i.e. if an ICB has a high admission rate for emergency elderly, they tend to have high admission rates in the frail cohort also.

```{r}
compare_2_2_plots("A. Frail", "B. Emergency elderly", "Frail", "Emergency elderly")
```

#### Frail x Falls

-   **Frail x Falls**: There is less evidence of correlation but almost all ICB's that have admission rates below 1 standard deviation in the falls cohorts also are below 1 standard deviation for frail admissions (lower left quadrant).

```{r}
compare_2_2_plots("A. Frail", "C. Falls", "Frail", "Falls")
```

#### Frail x End of life

-   **Frail x End of life**: There is no evidence of a link between an ICB's admission rate in end of life admissions and frail patients.

```{r}
compare_2_2_plots("A. Frail", "D. End of life", "Frail", "End of life")
```

#### Emergency elderly x Falls

-   **Emergency elderly x Falls**: Correlation is demonstrated again when considering falls and emergency elderly however the 4 ICB's with the highest admission rates for falls are within 1 standard deviation of the mean emergency elderly admission rate and the 5 of the highest 7th emergency elderly rates are similarly closer to average for the falls cohorts.  

```{r}
compare_2_2_plots("B. Emergency elderly", "C. Falls", "Emergency elderly", "Falls")
```

#### Emergency elderly x End of life

-   **Emergency elderly x End of life**: No evidence of correlation between end of life and emergency elderly admission rates with a number of ICB's position far from average on 1 but near the average on the other.

```{r}
compare_2_2_plots("B. Emergency elderly", "D. End of life", "Emergency elderly", "End of life")
```

#### Falls x End of life

-   **Falls x End of life**: Similarly, no evidence of correlation between end of life and falls admission rates with a number of ICB's position far from average on 1 but near the average on the other.

```{r}
compare_2_2_plots("C. Falls", "D. End of life", "Falls", "End of life")
```
:::

## Variation in outcomes

Furthering our assessment of variation in treatment of these mitigable patient populations, we analyse variation in selected patient outcomes; readmission and mortality.

This sub-strand of our descriptive analysis aims to investigate the risk of death or subsequent readmission in patients who have had an emergency hospital admission. We apply survival analysis techniques to help us understand what factors influence the risk of death or readmission, including determining variability between ICB area.

Full survival analysis output found here: <https://the-strategy-unit.github.io/Community_Strategies/community_strategy_survival_analysis.html>

![Cox regression example](cox_prop_elderly_readmission.png)

<hr style="border-top: 2px solid #f9bf07">

# Mitigator overlap

Finally, explore the correlation and overlap between and within the cohorts included within the Community Strategies analysis.

```{r}
compare_admission_rates |> 
  ggpairs(columns = 2:5, 
          upper = list(continuous = wrap("cor", size = 3)),
          lower = list(continuous = wrap("points", alpha = 0.6, size = 1.5)),
          diag = list(continuous = wrap("densityDiag", alpha = 0.5))) +
  theme_SU() +
  theme(strip.background = element_rect(fill = NA, colour = "grey")) +
  labs(title = "Admission Rates per ICB for mitigable sub-cohorts",
       subtitle = "SUS adjusted admission rate per 1,000 | 2023",
       x = "Admission rate")
```

Full analysis on overlap found here: <https://the-strategy-unit.github.io/Community_Strategies/identifying_overlap_between_cohorts.html>

![Overlap Venn diagram example](overlap_venn_plot.png)

# Contact

If you have any questions or comments regarding any of the above analysis please email [strategy.unit\@nhs.net](mailto:strategy.unit@nhs.net){.email} citing the NHSE Community Services Analysis.