Scotland disasters emergencies

R/capacity/disasters-emergencies/scotland/build-index.R

@@ -0,0 +1,181 @@
+# Fire station response time is not included for now (17/03) as don't have LTLA level data 
+# (only FRA level). Include again when get this from Gov. 
+# ---- Load libraries and Functions ----
+library(tidyverse)
+library(demographr)
+library(geographr)
+
+source("R/utils.R")
+
+# Load indicators ----
+indicators <-
+  load_indicators(
+    path = "data/capacity/disasters-emergencies/scotland",
+    key = "lad_code"
+  )
+
+# Check entry for all 2020 LAD codes
+lads_20 <- lookup_postcode_oa11_lsoa11_msoa11_ltla20 |>
+   filter(str_detect(ltla20_code, "^S")) |>
+   distinct(lad_code = ltla20_code) 
+
+lads_20 |>
+  anti_join(indicators, by = "lad_code")
+
+indicators |>
+  anti_join(lads_20, by = "lad_code")
+
+# Check for NAs
+indicators |>
+  dplyr::filter(if_any(everything(), ~is.na(.x))) 
+
+# Align direction so that high score = low capacity 
+indicators_aligned <- indicators |>
+  mutate(cap_exp_person = cap_exp_person * -1)
+
+# Check normality of indicators
+indicators_aligned |>
+  pivot_longer(civic_assets_extent:cap_exp_person, names_to = "variable", values_to = "value") %>%
+  ggplot(aes(x = value)) +
+  geom_density() +
+  facet_wrap(vars(variable), ncol = 3, scales = "free")
+
+for(i in 2:4) {
+  print(colnames(indicators[, i]))
+  print(shapiro.test(indicators[[i]])$p.value)
+  if (shapiro.test(indicators[[i]])$p.value < 0.05) {
+    print("Not normally distributed")
+  }
+  else {
+    print("Normally distributed")
+  }
+  cat("\n")  
+
+}
+
+#all the variables are not normally distributed
+
+# Exponential transformation ----
+
+ranked <- indicators_aligned |>
+  mutate_if(is.numeric, rank)
+
+scale01 <- function(x) (x - min(x))/diff(range(x))
+
+scale_ind <- ranked |>
+  mutate_if(is.numeric, scale01)
+
+exponential = function(x) (-23*log(1-x*(1-exp(1)^(-100/23))))
+
+exp_ind <- scale_ind |>
+  mutate_if(is.numeric, exponential)
+
+for(i in 2:4) {
+  print(colnames(exp_ind[, i]))
+  print(shapiro.test(exp_ind[[i]])$p.value)
+  if (shapiro.test(exp_ind[[i]])$p.value < 0.05) {
+    print("Not normally distributed")
+  }
+  else {
+    print("Normally distributed")
+  }
+  cat("\n")  
+
+}
+
+library(psych)
+for(i in 2:4) {
+  print(colnames(exp_ind[, i]))
+  print(skew(exp_ind[[i]]))
+  cat("\n")  
+
+}
+
+transformed <- exp_ind |>
+  mutate(civic_assets_extent = log10(civic_assets_extent+1),
+         engagement_extent = log10(engagement_extent+1),
+         cap_exp_person = log10(cap_exp_person+1))
+
+for(i in 2:4) {
+  print(colnames(transformed[, i]))
+  print(shapiro.test(transformed[[i]])$p.value)
+  if (shapiro.test(transformed[[i]])$p.value < 0.05) {
+    print("Not normally distributed")
+  }
+  else {
+    print("Normally distributed")
+  }
+  cat("\n")  
+
+}
+
+# After the transformation engagement_extent is still not normally distributed
+
+# Normalise indicators
+normalised <- transformed |>
+  normalise_indicators(ignore_nas = T)
+
+normalised |>
+  pivot_longer(civic_assets_extent:cap_exp_person, names_to = "variable", values_to = "value") %>%
+  ggplot(aes(x = value)) +
+  geom_density() +
+  facet_wrap(vars(variable), ncol = 3, scales = "free")
+
+# MFLA
+standardised = function(x) (x - mean(x))/sd(x)
+rank2 = function(x) rank(x, na.last = FALSE)
+
+mfla_score <- function(d) {
+
+  # Rank and normalise indicators to mean 0, SD 1.
+  d <- d %>%
+    mutate_if(is.numeric, list(scaled = function(x) standardised(rank2(x))))
+
+  # Extract weights
+  d_weights <- d %>%
+    select(ends_with("_scaled")) %>%
+    factanal(factors = 1) %>%
+    tidy() %>%
+    select(-uniqueness, weights = fl1) %>%
+    mutate(weights = abs(weights),
+           weights = weights/sum(weights))
+
+  # Multiply model weights by respective column to get weighted indicators
+  d_weighted_ind <- d %>%
+    select(d_weights$variable) %>%
+    map2_dfc(d_weights$weights, `*`) %>%
+    select_all(list(~ str_remove(., "_scaled"))) %>%
+    select_all(list(~ str_c(., "_weighted")))
+
+  # Combine weighted indicators with original data
+  d <- bind_cols(d, d_weighted_ind)
+
+  # Sum weighted indicators
+  d <- d %>%
+    mutate(mfla_score = reduce(select(., ends_with("_weighted")), `+`))
+
+  # Return data
+  return(d)
+
+}
+
+#Apply the function
+indicators_mfla <- normalised |>
+  mfla_score() |>
+  select(lad_code,
+         mfla_score)
+
+# Domain scores
+indicators_scores <- indicators_mfla |>
+  calculate_domain_scores(domain_name = "capacity")
+
+# Inverting ranks and deciles so that higher scores = higher capacity
+indicators_invert <- indicators_scores |>
+  mutate(capacity_domain_rank = inverse_rank(capacity_domain_rank),
+         capacity_domain_quantiles = invert_this(capacity_domain_quantiles)) |>
+  select(lad_code,
+         deciles = capacity_domain_quantiles)
+
+# Save index
+indicators_invert |>
+  write_csv("data/capacity/disasters-emergencies/scotland/capacity-index.csv")

R/capacity/disasters-emergencies/scotland/community-assets.R

@@ -0,0 +1,95 @@
+# Load packages
+library(tidyverse)
+library(readxl)
+library(geographr)
+library(demographr)
+
+source("R/utils.R")
+
+# Community Needs Index----
+raw <- read_excel("data/on-disk/SCNI and domains - 2021.xlsx")
+
+civic_assets <- raw |>
+  select(
+    msoa_code = "Area code",
+    msoa_name = "Area name",
+    civic_assets_score = "Civic_Assets_Domain",
+    civic_assets_rank = "Civic_Assets_Domain_rank"
+  ) |>
+  mutate(msoa_code = gsub('s','S', msoa_code))
+
+# Population data -----
+pop_msoa <- population_dz_20 |>
+  filter_codes(dz_code, "^S") |>
+  filter(sex == "All") |>
+  select(dz_code,
+         total_population) |>
+  inner_join(lookup_postcode_oa11_lsoa11_msoa11_ltla20 |>
+               filter(str_detect(msoa11_code, "^S")) |>
+               select(msoa_code = msoa11_code,
+                      dz_code = lsoa11_code) |>
+               distinct()) |>
+  group_by(msoa_code) |>
+  summarise(total_population = mean(total_population))
+
+# Check MSOA with no population data
+civic_assets |>
+  left_join(pop_msoa, by = "msoa_code") |>
+  filter(is.na(total_population))
+
+# Lookup codes ----
+lookup_lad <- lookup_postcode_oa11_lsoa11_msoa11_ltla20 |>
+  filter(str_detect(msoa11_code, "^S")) |>
+  select(msoa_code = msoa11_code,
+         lad_code = ltla20_code) |>
+  distinct()
+
+# Check any missing LADs
+lookup_lad |>
+  anti_join(civic_assets, by = "msoa_code")
+# 0
+civic_assets |>
+  anti_join(lookup_lad, by = "msoa_code")
+#0
+
+# Calculate extent----
+# HIGH SCORE/LOW RANK = LOW CAPABILITIES
+civic_assets_lad <- civic_assets |>
+  left_join(lookup_lad, by = "msoa_code") |>
+  left_join(pop_msoa, by = "msoa_code") |>
+  calculate_extent(
+    var = civic_assets_score,
+    higher_level_geography = lad_code,
+    population = total_population,
+    weight_high_scores = TRUE #  TRUE when a highest variable score equates to a lower capacity
+  ) |>
+  rename(civic_assets_extent = extent)
+
+# Lookup names (for graph readability)
+lookup_names <- lookup_ltla_ltla |>
+  filter(str_detect(ltla20_code, "^S")) |>
+  select(lad_name = ltla20_name,
+         lad_code = ltla20_code) |>
+  distinct()
+
+civic_assets_lad |>
+  left_join(lookup_names, by ="lad_code") |>
+  mutate(lad_name = fct_reorder(lad_name, desc(civic_assets_extent))) |>
+  ggplot(aes(x = lad_name, y = civic_assets_extent))+
+  geom_point() +
+  theme_classic() +
+  labs(title = "Civic Assets Extent by LAD",
+       x = "LAD",
+       y = "Civic Assets Extent") +
+  guides(x = guide_axis(angle = 90))
+
+civic_assets_lad |>
+  ggplot(aes(y = civic_assets_extent)) +
+  geom_boxplot() +
+  ylab("Civic Assets Extent") +
+  theme(axis.text.x = element_blank(),
+        axis.ticks.x = element_blank())
+
+# Save data -----
+civic_assets_lad |>
+  write_rds("data/capacity/disasters-emergencies/scotland/community-assets.rds")

R/capacity/disasters-emergencies/scotland/engagement.R

@@ -0,0 +1,99 @@
+# Load packages
+library(tidyverse)
+library(readxl)
+library(geographr)
+library(demographr)
+
+source("R/utils.R")
+
+
+# OSCI Community Needs Index data ----
+raw <- read_excel("data/on-disk/SCNI and domains - 2021.xlsx")
+
+engagement <- raw |>
+  select(
+    msoa_code = "Area code",
+    msoa_name = "Area name",
+    engagement_score = "ActiveEngaged_Domain",
+    engagement_rank = "ActiveEngaged_Domain_rank"
+  ) |>
+  mutate(msoa_code = gsub('s','S', msoa_code))
+
+cor(engagement$engagement_score, engagement$engagement_rank)
+
+#  Population data -----
+pop_msoa <- population_dz_20 |>
+  filter_codes(dz_code, "^S") |>
+  filter(sex == "All") |>
+  select(dz_code,
+         total_population) |>
+  inner_join(lookup_postcode_oa11_lsoa11_msoa11_ltla20 |>
+               filter(str_detect(msoa11_code, "^S")) |>
+               select(msoa_code = msoa11_code,
+                      dz_code = lsoa11_code) |>
+               distinct()) |>
+  group_by(msoa_code) |>
+  summarise(total_population = mean(total_population))
+
+# Check MSOA with no population data
+engagement |>
+  left_join(pop_msoa, by = "msoa_code") |>
+  filter(is.na(total_population))
+
+# Lookup codes ----
+lookup_lad <- lookup_postcode_oa11_lsoa11_msoa11_ltla20 |>
+  filter(str_detect(msoa11_code, "^S")) |>
+  select(msoa_code = msoa11_code,
+         lad_code = ltla20_code) |>
+  distinct()
+
+# Check any missing LADs
+lookup_lad |>
+  anti_join(engagement, by = "msoa_code")
+# 0
+engagement |>
+  anti_join(lookup_lad, by = "msoa_code")
+#0
+
+# Calculate extent----
+# HIGH SCORE/LOW RANK = LOW CAPABILITIES
+
+engagement_lad <- engagement |>
+  left_join(lookup_lad, by = "msoa_code") |>
+  left_join(pop_msoa, by = "msoa_code") |>
+  calculate_extent(
+    var = engagement_score,
+    higher_level_geography = lad_code,
+    population = total_population,
+    weight_high_scores = TRUE #  TRUE when a highest variable score equates to a lower capacity
+  ) |>
+  rename(engagement_extent = extent)
+
+# Lookup names (for graph readability)
+lookup_names <- lookup_ltla_ltla |>
+  filter(str_detect(ltla20_code, "^S")) |>
+  select(lad_name = ltla20_name,
+         lad_code = ltla20_code) |>
+  distinct()
+
+engagement_lad |>
+  left_join(lookup_names, by ="lad_code") |>
+  mutate(lad_name = fct_reorder(lad_name, desc(engagement_extent))) |>
+  ggplot(aes(x = lad_name, y = engagement_extent))+
+  geom_point() +
+  theme_classic() +
+  labs(title = "Engagement Extent by LAD",
+       x = "LAD",
+       y = "Engagement Extent") +
+  guides(x = guide_axis(angle = 90))
+
+engagement_lad |>
+  ggplot(aes(y = engagement_extent)) +
+  geom_boxplot() +
+  ylab("Engagement Extent") +
+  theme(axis.text.x = element_blank(),
+        axis.ticks.x = element_blank())
+
+# Save data -----
+engagement_lad |>
+  write_rds("data/capacity/disasters-emergencies/scotland/engagement.rds")