Last active
December 5, 2022 16:46
-
-
Save rcatlord/a09e118a8194993cd0d0b219735acd75 to your computer and use it in GitHub Desktop.
Measuring spatial autocorrelation
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# 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') |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment