create vulnerability index - LA.r

##
## Aggregate England's MSOA-level Vulnerability Index into Local Authorities
##
library(tidyverse)
library(janitor)
library(sf)

source("functions.r")
source("load lookup tables.r")

# run "create vulnerability index - MSOA.r" to make this dataset
vi = read_csv("output/vulnerability-MSOA-UK.csv")


##
## lookup tables
##
# MSOA to LA codes for England, Wales and Scotland
msoa_lad = load_lookup_lsoa_msoa_lad() %>% 
  select(MSOA11CD, LAD17CD) %>% 
  distinct()

# SOA to LA codes for Northern Ireland
soa_lgd = load_lookup_sa_lgd() %>% 
  select(MSOA11CD = LSOA11CD, LAD17CD = LAD18CD) %>% 
  distinct()

# append NI's codes to the rest
msoa_lad = bind_rows(msoa_lad, soa_lgd)

lad_17_19 = read_csv("data/LAD 2017 to LAD 2019 codes.csv")  # lookup of Local Authority codes from 2017 to 2019

lad_names = read_csv("https://opendata.arcgis.com/datasets/c3ddcd23a15c4d7985d8b36f1344b1db_0.csv") %>% 
  select(LAD19CD, LAD19NM)

# Local Authority District to Local Resilience Forum (December 2019) Lookup in England and Wales
# https://geoportal.statistics.gov.uk/datasets/local-authority-district-to-local-resilience-forum-december-2019-lookup-in-england-and-wales
# lad_lrf = read_csv("https://opendata.arcgis.com/datasets/203ea94d324b4fcda24875e83e577060_0.csv")

# MSOA to LRF lookup
# msoa_lrf = msoa_lad %>% 
#   left_join(lad_17_19, by = "LAD17CD") %>% 
#   left_join(lad_lrf, by = "LAD19CD") %>% 
#   select(MSOA11CD, LRF19CD) %>% 
#   distinct()

##
## get MSOA-level population estimates
##
source("prep population estimates.r")


##
## aggregate vulnerability domains into Local Authorities by calculating:
## - proportions of 20% most-vulnerable neighbourhoods in each LA
## - extent of the population living in the most-vulnerable neighbourhoods
## - population-weighted average vulnerability score
##
vi_lad = vi %>% 
  left_join(msoa_lad, by = c("Code" = "MSOA11CD")) %>% 
  left_join(lad_17_19, by = "LAD17CD") %>% 
  left_join(pop, by = "Code")

# clinical vulnerability
clinical = vi_lad %>% 
  aggregate_scores(domain = "Clinical", aggregate_by = "LAD19CD") %>% 
  rename(`Proportion of neighbourhoods in 20% most clinically vulnerable` = Proportion,
         `Extent of population living in highly clinically vulnerable areas` = Extent,
         `Population-weighted clinical vulnerability score` = Score)

# health/wellbeing vulnerability 
health_wellbeing = vi_lad %>% 
  aggregate_scores(domain = "Health/Wellbeing", aggregate_by = "LAD19CD") %>% 
  rename(`Proportion of neighbourhoods in 20% most health/wellbeing vulnerable` = Proportion,
         `Extent of population living in highly health/wellbeing vulnerable areas` = Extent,
         `Population-weighted health/wellbeing vulnerability score` = Score)

# economic vulnerability
economic = vi_lad %>% 
  aggregate_scores(domain = "Economic", aggregate_by = "LAD19CD") %>% 
  rename(`Proportion of neighbourhoods in 20% most economically vulnerable` = Proportion,
         `Extent of population living in highly economically vulnerable areas` = Extent,
         `Population-weighted economic vulnerability score` = Score)

# social vulnerability 
social = vi_lad %>% 
  aggregate_scores(domain = "Social", aggregate_by = "LAD19CD") %>% 
  rename(`Proportion of neighbourhoods in 20% most socially vulnerable` = Proportion,
         `Extent of population living in highly socially vulnerable areas` = Extent,
         `Population-weighted social vulnerability score` = Score)

# socioeconomic vulnerability 
socioeconomic = vi_lad %>% 
  aggregate_scores(domain = "Socioeconomic", aggregate_by = "LAD19CD") %>% 
  rename(`Proportion of neighbourhoods in 20% most socioeconomically vulnerable` = Proportion,
         `Extent of population living in highly socioeconomically vulnerable areas` = Extent,
         `Population-weighted socioeconomic vulnerability score` = Score)

# overall vulnerability (without the LA-level indicators added below)
overall = vi_lad %>% 
  aggregate_scores(aggregate_by = "LAD19CD") %>% 
  rename(`Proportion of neighbourhoods in 20% most vulnerable` = Proportion,
         `Extent of population living in highly vulnerable areas` = Extent,
         `Population-weighted vulnerability score` = Score)

##
## Include LA-level indicators
##
asylum = read_csv("data/asylum-LA.csv") %>% 
  select(LAD19CD, `People receiving Section 95 support`)
  # mutate(`People receiving Section 95 support` = replace_na(`People receiving Section 95 support`, 0))

community_needs_lad = read_csv("data/community-needs-LA.csv") %>%
  select(-`Proportion of wards with greatest community needs`)


##
## Calculate new overall vulnerability index
##
vulnerability = clinical %>%
  select(LAD19CD, `Population-weighted clinical vulnerability score`) %>% 
  
  left_join(health_wellbeing %>% select(LAD19CD, starts_with("Population-weighted")), 
            by = "LAD19CD") %>% 
  
  left_join(economic %>% select(LAD19CD, starts_with("Population-weighted")), 
                                by = "LAD19CD") %>% 
  
  left_join(social %>% select(LAD19CD, ends_with("Vulnerability score")),
            by = "LAD19CD") %>% 
  
  left_join(asylum, by = "LAD19CD") %>% 
  mutate(`People receiving Section 95 support` = replace_na(`People receiving Section 95 support`, 0)) %>%  # missing values mean no people seeking asylum in that LA
  
  left_join(community_needs_lad, by = "LAD19CD") %>% 
  
  mutate(Country = get_country(LAD19CD))

# helper function for calculating within-nation vulnerability scores
calc_LA_scores = function(d) {
  d %>% 
    # calculate ranks
    mutate(
      `Clinical Vulnerability rank` = rank2(`Population-weighted clinical vulnerability score`),
      `Health/Wellbeing Vulnerability rank` = rank2(`Population-weighted health/wellbeing vulnerability score`),
      `Economic Vulnerability rank` = rank2(`Population-weighted economic vulnerability score`),
      `Social Vulnerability rank` = rank2(`Population-weighted social vulnerability score`),
      
      `Aslym rank` = rank2(standardised(rank2(`People receiving Section 95 support`))),
      `Community Assets rank` = rank2(standardised(rank2(`Proportion of wards with worst Community Assets scores`))),
      `Community Engagement rank` = rank2(standardised(rank2(`Proportion of wards with worst Engagement scores`))),
    ) %>% 
    select(LAD19CD, Country, ends_with("rank")) %>% 
    
    calc_domain_scores(rank.indicators = FALSE, bespoke.domains = TRUE)
}

# split into separate datasets for each country then calculate vulnerability scores - until can work out how to get this code working:
#   group_by(Country) %>% 
#   calc_domain_scores() %>% 
#   ungroup() %>% 
v_eng = vulnerability %>% 
  filter(Country == "England") %>% 
  calc_LA_scores()

v_wal = vulnerability %>% 
  filter(Country == "Wales") %>% 
  calc_LA_scores()

v_sco = vulnerability %>% 
  filter(Country == "Scotland") %>% 
  calc_LA_scores()

v_ni = vulnerability %>% 
  filter(Country == "Northern Ireland") %>% 
  calc_LA_scores()

# recombine
vulnerability = bind_rows(v_eng, v_wal, v_sco, v_ni) %>% 
  select(LAD19CD, starts_with("Vulnerability"))

# merge all proportions togethers
vulnerability_all = clinical %>% 
  left_join(health_wellbeing, by = "LAD19CD") %>% 
  left_join(economic, by = "LAD19CD") %>% 
  left_join(social, by = "LAD19CD") %>% 
  left_join(socioeconomic, by = "LAD19CD") %>% 
  left_join(overall, by = "LAD19CD") %>% 
  
  left_join(asylum, by = "LAD19CD") %>% 
  mutate(`People receiving Section 95 support` = replace_na(`People receiving Section 95 support`, 0)) %>%  # missing values mean no people seeking asylum in that LA
  left_join(community_needs_lad, by = "LAD19CD") %>% 
  
  left_join(vulnerability, by = "LAD19CD") %>% 
  
  mutate(Country = get_country(LAD19CD)) %>% 
  
  # calculate quintiles
  group_by(Country) %>% 
  mutate(`Clinical Vulnerability quintile` = calc_risk_quantiles(`Population-weighted clinical vulnerability score`, quants = 5),
         `Health/Wellbeing Vulnerability quintile` = calc_risk_quantiles(`Population-weighted health/wellbeing vulnerability score`, quants = 5),
         `Economic Vulnerability quintile` = calc_risk_quantiles(`Population-weighted economic vulnerability score`, quants = 5),
         `Social Vulnerability quintile` = calc_risk_quantiles(`Population-weighted social vulnerability score`, quants = 5),
         `Socioeconomic Vulnerability quintile` = calc_risk_quantiles(`Population-weighted socioeconomic vulnerability score`, quants = 5),
         `Vulnerability quintile` = calc_risk_quantiles(`Vulnerability rank`, quants = 5)) %>% 
  ungroup() %>% 
  
  # get LA names
  left_join(lad_names, by = "LAD19CD") %>% 
  select(Code = LAD19CD, Name = LAD19NM, everything())

# calculate quintiles
vulnerability_scores = vulnerability_all %>% 
  select(Code, Name, ends_with("quintile"))


##
## save
##
vulnerability_all %>% write_csv("output/vulnerability-LA.csv")

vulnerability_all %>% 
  select(Code, Name, ends_with("quintile")) %>% 
  write_csv("output/vulnerability-scores-LA.csv")

# Local Authority Districts (December 2019) Boundaries UK BUC
# source: https://geoportal.statistics.gov.uk/datasets/local-authority-districts-december-2019-boundaries-uk-buc
lads = read_sf("https://opendata.arcgis.com/datasets/3a4fa2ce68f642e399b4de07643eeed3_0.geojson")

geojson_name = "output/vulnerability-LA.geojson"
if (file.exists(geojson_name)) file.remove(geojson_name)

lads %>% 
  left_join(vulnerability_all, by = c("lad19cd" = "Code")) %>% 
  
  # append a note to deciles and ranks
  rename_at(vars(matches("decile|rank")), ~ str_c(., " (1 = least vulnerable)")) %>% 
  
  select(Code = lad19cd, Name, everything(), -objectid, -lad19nm, -lad19nmw, -st_areashape, -st_lengthshape, -bng_e, -bng_n, -long, -lat) %>% 
  write_sf(geojson_name)