Measuring spatial autocorrelation

README.md

script.R

# Measuring spatial autocorrelation #

library(tidyverse) ; library(httr) ; library(readxl) ; library(sf) ; library(spdep) ; library(magrittr) ; library(tmap) ; library(htmltools) ; library(leaflet)

# Load data ------------------------------------------

# Income deprivation
# Source: ONS
# URL: https://www.ons.gov.uk/peoplepopulationandcommunity/personalandhouseholdfinances/incomeandwealth/datasets/mappingincomedeprivationatalocalauthoritylevel
tmp <- tempfile(fileext = ".xlsx")
GET(url = "https://www.ons.gov.uk/file?uri=/peoplepopulationandcommunity/personalandhouseholdfinances/incomeandwealth/datasets/mappingincomedeprivationatalocalauthoritylevel/2019/localincomedeprivationdata.xlsx",
    write_disk(tmp))
df <- read_xlsx(tmp, sheet = "LSOA") %>% 
  select(lsoa11cd = `LSOA code (2011)`,
         lsoa11nm = `LSOA name (2011)`,
         ltla19cd = `Local Authority District code (2019)`,
         ltla19nm = `Local Authority District name (2019)`,
         score = `Income Score (rate)`,
         rank = `Income Rank (where 1 is most deprived)`)

# LSOA boundaries
# Source: ONS Open Geography Portal
# URL: https://geoportal.statistics.gov.uk/datasets/ons::lsoa-dec-2011-boundaries-super-generalised-clipped-bsc-ew-v3/explore
sf <- st_read("https://services1.arcgis.com/ESMARspQHYMw9BZ9/arcgis/rest/services/LSOA_Dec_2011_Boundaries_Super_Generalised_Clipped_BSC_EW_V3_2022/FeatureServer/0/query?outFields=*&where=1%3D1&f=geojson") %>% 
  select(lsoa11cd = LSOA11CD)

# select a local authority ------------------------------------------
id <- "Liverpool"

# join, filter and covert data ------------------------------------------
sp <- left_join(sf, df, by = "lsoa11cd") %>% 
  filter(ltla19nm == id) %>%
  as("Spatial")

# Global Moran's I ------------------------------------------

# create neighbours list (first order queen contiguity)
nb <- poly2nb(sp, queen = TRUE, row.names = rownames(sp@data))

# calculate Global Moran's I
global_moran <- moran.test(sp$score, listw = nb2listw(nb, style = "W"))
print(paste(id, "has a Global Moran's Score of", 
            round(global_moran[["estimate"]][["Moran I statistic"]],2)))

# create Moran's I scatterplot
moran.plot(sp$score, listw = nb2listw(nb, style = "W"), labels = sp$lsoa11nm,
           xlab = "Income deprivation score", ylab = "Lagged Income deprivation score")

# Local Moran's I ------------------------------------------

# calculate LISA scores
local_moran <- localmoran(x = sp$score, listw = nb2listw(nb, style = "W"), zero.policy = TRUE, na.action = na.omit)
significanceLevel <- 0.01
meanVal <- mean(sp$score)
local_moran %<>% as_tibble() %>%
  set_colnames(c("Ii","E.Ii","Var.Ii","Z.Ii","Pr(z > 0)")) %>%
  mutate(quad_sig = case_when(
    `Pr(z > 0)` > 0.05 ~ "Not significant",
    `Pr(z > 0)` <= 0.05 & Ii >= 0 & sp$score >= meanVal ~ "High-High",
    `Pr(z > 0)` <= 0.05 & Ii >= 0 & sp$score < meanVal ~ "Low-Low",
    `Pr(z > 0)` <= 0.05 & Ii < 0 & sp$score >= meanVal ~ "High-Low",
    `Pr(z > 0)` <= 0.05 & Ii < 0 & sp$score < meanVal ~ "Low-High")
  )
sp$quad_sig <- local_moran$quad_sig %>% replace_na("Not significant")

lisa <- st_as_sf(sp) %>% 
  mutate(popup = str_c("<strong>", lsoa11nm, "</strong><br/>", quad_sig, "<br/>Income score: ", score, "<br/>Income rank: ", format(rank, big.mark = ",")) %>% map(HTML))

# static map
tmap_mode("plot") # or view
tm <- tm_shape(lisa) +
  tm_polygons("quad_sig", 
              palette = c("Not significant" = "#F0F0F0", 
                          "High-High" = "#E93F36", 
                          "Low-Low" = "#2144F5", 
                          "Low-High" = "#9794F8", 
                          "High-Low" = "#EF9493"),
              title = "LISA Cluster Map",
              border.col = "#000000",
              legend.hist = TRUE) + 
  tm_layout(frame = FALSE,
            main.title = paste("Income deprivation in", id),
            main.title.position = "left",
            main.title.size = 1.2,
            main.title.color = "#212121",
            legend.outside = TRUE,
            legend.hist.width = 5) +
  tm_credits(paste("Global Moran's I:", round(global_moran[["estimate"]][["Moran I statistic"]],2)), 
             fontface = "bold",
             align = "left",
             position = c("left", "bottom"))
tm
# tmap_save(tm, "LISA_map.png", height = 7)

# interactive map
factpal <- colorFactor(c("#F0F0F0", "#E93F36", "#2144F5", "#9794F8", "#EF9493"),
                       levels = c("Not significant", "High-High", "Low-Low", "Low-High", "High-Low"),
                       ordered = TRUE)
leaflet() %>% 
  addProviderTiles("CartoDB.Positron") %>% 
  addPolygons(data = lisa, fillColor = ~factpal(quad_sig), fillOpacity = 0.6, 
              stroke = TRUE, color = "#212121", weight = 1, layerId = ~lsoa11cd,
              label = ~popup, 
              highlight = highlightOptions(color = "#FFFF00", weight = 3, opacity = 1, bringToFront = TRUE)) %>%
  addLegend(position = "bottomleft", colors = c("#f0f0f0", "#E93F36", "#2144F5", "#9794F8", "#EF9493"), opacity = 0.6,
            labels = c("Non-sig", "High-High", "Low-Low", "Low-High", "High-Low"),
            title = "LISA Cluster Map") %>% 
  addControl(paste0("<strong>Income deprivation in ", id, "</strong>"), position = 'topright')