How can I create a LISA cluster map in R?

[LucasGBegiato]

I've read in many places how to create a LISA map, but I'm not really understanding the process. I already have the SHAPEFILE and the DATA SET together, I would like to know how do I get a figure of the following type using the "Incidência da COVID-19" variable resulted after I "full_joined" to variable "Data".

--------- Packages -------------------------------

setwd("C:/Users/lucas/Desktop/Relatório Final FAPESP")

if (!require("geobr")) install.packages("geobr")
if (!require("tidyverse")) install.packages("tidyverse")
if (!require("ggspatial")) install.packages("ggspatial")
if (!require("magrittr")) install.packages("magrittr")
if (!require("lubridate")) install.packages("lubridate")
if (!require("readxl")) install.packages("readxl")
library(readxl)
library(geobr)
library(tidyverse)
library(ggspatial)
library(magrittr)
library(lubridate)

--------- Shape File ------------------------------

SP.Municipios <- read_municipality(code_muni = "all",year=2019) %>%
filter(code_state == 35) %>%
select(-c("code_state","abbrev_state","code_region","name_region","name_state")) %>%
rename(Municípios = name_muni)

ggplot(SP.Municipios) +
geom_sf()

--------- Data Set -----------------------------

Data.Raw <- read_excel("Base.xlsx")

Total.Data <- full_join(Data.Raw, SP.Municipios, by = "Municípios")

Data <- Total.Data[c(-644,-645),]

------------------------------------

the DATA is like this:

Can anyone help me with this? I don't think it's that difficult, but unfortunately I don't understand what I have to do. If anyone needs the database, let me know.

[jlacko]

You may want to read a bit on the subject; the process behind LISA (spatial lag, weights, Moran and local Moran) is heavy on theory.

But if your aim is to prepare a plot first, and improve on your understanding later - which is a perfectly legit aim! - you can consider something along these lines of code.

It relies heavily on {sfweight} package by Josh Parry, who might be active on this site. It provides a friendlier interface to older, and somewhat more general, but way more complex, {spdep} package. I am not certain if it is on CRAN at present, if not you can use the GitHub version.

As I don't have access to your data I am using the trusty old North Carolina shapefile that ships with {sf} package; you can easily replace it with your dataset.

You need to have three things in place before analyzing LISA clusters:

• spatial neighbors (read up on queen and rook approach; queen is the default); implemented via sfweight::st_neighbors()
• spatial weights, implemented via sfweight::st_weights() with the neighbors as argument
• spatial lag of a variable, implemented via sfweight::st_lag(), with neigbors and weights as additional arguments

With these in place you can calculate the LISA categories using sfweight::categorize_lisa(), the arguments are your variable of interest and its spatial lag.

As I wrote, it is not the full story, but it should give you a start now to improve upon later...

library(sf)
library(dplyr)
library(sfweight) # remotes::install_github("Josiahparry/sfweight")
library(ggplot2)

shape <- st_read(system.file("shape/nc.shp", package="sf")) # included with sf package

# calucualte the lisa groups
shape_lisa <- shape %>% 
  mutate(nb = st_neighbors(geometry),
         wts = st_weights(nb),
         lag_SID79 = st_lag(SID79, nb, wts),
         lisa = categorize_lisa(SID79, lag_SID79))

# report results
ggplot(data = shape_lisa) +
  geom_sf(aes(fill = lisa))

[LucasGBegiato]

Hi, you again! Thank you for helping me. I understand everything you explained, including the code! But mine was wrong because of an empty neighborhood (just like what happened in the SAR Model). So, thanks to your help back there, I was able to work around this problem and fix the empty neighborhood error!

--------- Packages ---------

setwd("C:/Users/lucas/Desktop/Relatório Final FAPESP")

if (!require("geobr")) install.packages("geobr")
if (!require("tidyverse")) install.packages("tidyverse")
if (!require("ggspatial")) install.packages("ggspatial")
if (!require("magrittr")) install.packages("magrittr")
if (!require("lubridate")) install.packages("lubridate")
if (!require("readxl")) install.packages("readxl")
library(readxl)
library(geobr)
library(tidyverse)
library(ggspatial)
library(magrittr)
library(lubridate)

--------- Shape File ---------

SP.Municipios <- read_municipality(code_muni = "all",year=2019) %>%
filter(code_state == 35) %>%
select(-c("code_state","abbrev_state","code_region","name_region","name_state")) %>%
rename(Municípios = name_muni) %>%

ggplot(SP.Municipios) +
geom_sf()

--------- Data Set ---------

Data.Raw <- read_excel("Base.xlsx")

Total.Data <- full_join(Data.Raw, SP.Municipios, by = "Municípios")

Data <- Total.Data[c(-644,-645),] %>%
rename(Incid = Incidência da COVID-19)

Data <- st_as_sf(Data)

--------- Clusters Analysis ---------

if (!require("remotes")) install.packages("remotes")
library(remotes)

remotes::install_github("Josiahparry/sfweight")

if (!require("sf")) install.packages("sf")
if (!require("dplyr")) install.packages("dplyr")
if (!require("ggplot2")) install.packages("ggplot2")
library(sf)
library(dplyr)
library(ggplot2)
library(sfweight)
if (!require("readr")) install.packages("readr")
if (!require("rgdal")) install.packages("rgdal")
if (!require("spdep")) install.packages("spdep")
if (!require("spatialreg")) install.packages("spatialreg")
library(readr)
library(rgdal)
library(spdep)
library(spatialreg)

#This will cause an error.
wrong <- Data %>%
poly2nb() %>%
nb2listw()

problematic <- Data %>%
poly2nb() %>%
card() %>%
magrittr::equals(0) %>%
which()

corrected_Data <- Data[-problematic, ]

shape_lisa.SP <- corrected_Data %>%
mutate(nb = st_neighbors(geom),
wts = st_weights(nb),
lag_Incid = st_lag(Incid, nb, wts),
lisa = categorize_lisa(Incid, lag_Incid))

ggplot(data = shape_lisa.SP) +
geom_sf(aes(fill = lisa))

--------- / --------- / --------- / --------- / ---------

Would I be able to modify the caption, from that, adding the "not significant" to p<0.05? Because this map is a little visually polluted. Something just like that:

If it's not possible, that's totally fine. You helped me a lot and I am again extremely grateful for that.

Virtual hugs from Brazil,
Lucas.

[jlacko]

Glad to be of service; spatial analysis is fun!

The Pakistan labels kind of remind me of the results produced by GEODA; you can use it via the {rgeoda} package, directly from the guys in Chicago / Luc Anselin's lab.

Consider this code:

library(sf)
library(rgeoda)
library(ggplot2)

shape <- st_read(system.file("shape/nc.shp", package="sf")) # included with sf package

# create weights object
queen_w <- queen_weights(shape)

# calculate LISA as per GEODA
lisa <- local_moran(queen_w, shape["SID79"]) # or any other variable :)

# process results
shape$cluster <- as.factor(lisa$GetClusterIndicators())
levels(shape$cluster) <- lisa$GetLabels()

# A visual overview
ggplot(data = shape) +
  geom_sf(aes(fill = cluster))

You would want to tune the fill palette a bit, but that is standard ggplot stuff...

I raise an imaginary beer glass from Prague, in a general south westerly direction

J.

[LucasGBegiato]

You are awesome! I tested it with the shape you used and it worked, but for some reason I don't know, when I put my shape, it says it's not in the 'sf' type.

I don't really know what's happening.

I raise an imaginary water glass, in the north east direction!
Thank's again.
L.

[LucasGBegiato]

queen_w <- queen_weights(SP.Municipios)

lisa <- local_moran(queen_w, corrected_Data["Incid"])

corrected_Data$cluster <- as.factor(lisa$GetClusterIndicators())

levels(corrected_Data$cluster) <- lisa$GetLabels()

ggplot(data = corrected_Data) +
geom_sf(aes(fill = cluster))

#---------------------------------------------

See that. I think i got it. I used my original shape file, and after that, I included the corrected_Data in local_moran.

So, my next step is to change the palette of colors, change the names, ...
I'm going to research about changing this color palette.

Thanks again and you can disconsider my last reply.

[LucasGBegiato]

GOT IT!!!!!