Local Moran’s I and spatial scales (using R and tmap)

Hi there, in this brief tutorial I will show you how to create a lisa map from scratch in R… For some reason the p-values in R and the p-values computed in Geoda for the local Moran’s I are slightly different. Any questions (or suggestions in order to improve this tutorial) are always welcome! Regards, Felipe.

1 loading required packages

knitr::opts_chunk$set(warning = FALSE, message = FALSE)
library(tidyverse)

## -- Attaching packages --------

## v ggplot2 3.3.0     v purrr   0.3.4
## v tibble  3.0.1     v dplyr   0.8.5
## v tidyr   1.0.3     v stringr 1.4.0
## v readr   1.3.1     v forcats 0.5.0

## -- Conflicts -----------------
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()

library(sp)
library(sf)

## Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1

library(fastmap)
library(skimr)
library(tmap)
library(tmaptools)
library(spdep)

## Loading required package: spData

## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`

library(rgeos)

## rgeos version: 0.5-3, (SVN revision 634)
##  GEOS runtime version: 3.8.0-CAPI-1.13.1 
##  Linking to sp version: 1.4-1 
##  Polygon checking: TRUE

#remotes::install_github("lixun910/rgeoda")
library(rgeoda)

## 
## Attaching package: 'rgeoda'

## The following object is masked from 'package:spdep':
## 
##     skater

# Change the presentation of decimal numbers to 4 and avoid scientific notation
options(prompt="R> ", digits=7, scipen=999)

2 Import the different shapefiles for Pakistan

nat<- st_read("shapefiles/gadm36_PAK_0.shp")

## Reading layer `gadm36_PAK_0' from data source `C:\Users\felip\OneDrive\Desktop\Meidai\Seminar Research\GitHub\tutorial-local-morans-I-spatial-scales\shapefiles\gadm36_PAK_0.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 2 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 60.89944 ymin: 23.70292 xmax: 77.84308 ymax: 37.09701
## geographic CRS: WGS 84

nat.sp <- as(nat, "Spatial")
class(nat.sp)

## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"

nat2<- st_read("shapefiles/gadm36_PAK_1.shp")

## Reading layer `gadm36_PAK_1' from data source `C:\Users\felip\OneDrive\Desktop\Meidai\Seminar Research\GitHub\tutorial-local-morans-I-spatial-scales\shapefiles\gadm36_PAK_1.shp' using driver `ESRI Shapefile'
## Simple feature collection with 8 features and 10 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 60.89944 ymin: 23.70292 xmax: 77.84308 ymax: 37.09701
## geographic CRS: WGS 84

nat2.sp <- as(nat2, "Spatial")
class(nat2.sp)

## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"

pro<- st_read("shapefiles/gadm36_PAK_2noNA.shp")

## Reading layer `gadm36_PAK_2noNA' from data source `C:\Users\felip\OneDrive\Desktop\Meidai\Seminar Research\GitHub\tutorial-local-morans-I-spatial-scales\shapefiles\gadm36_PAK_2noNA.shp' using driver `ESRI Shapefile'
## Simple feature collection with 29 features and 14 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 60.89944 ymin: 23.70292 xmax: 75.367 ymax: 36.8898
## geographic CRS: WGS 84

pro.sp <- as(pro, "Spatial")
class(pro.sp)

## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"

dis<- st_read("shapefiles/gadm36_PAK_3_noNA.shp")

## Reading layer `gadm36_PAK_3_noNA' from data source `C:\Users\felip\OneDrive\Desktop\Meidai\Seminar Research\GitHub\tutorial-local-morans-I-spatial-scales\shapefiles\gadm36_PAK_3_noNA.shp' using driver `ESRI Shapefile'
## Simple feature collection with 116 features and 21 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 60.89944 ymin: 23.70292 xmax: 75.367 ymax: 36.8898
## geographic CRS: WGS 84

#dis<- st_read("shapefiles/gadm36_PAK_3.shp")
#dis<- dis %>% filter(!is.na(ANC))
dis.sp <- as(dis, "Spatial")
class(dis.sp)

## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"

quick-maps

qtm(nat)

qtm(pro)

qtm(dis)

mun_merge<- dis.sp
#st_drop_geometry(dis)
#mun_merge$ANC

2.1 Do the shapefiles have the same projection?

proj4string(mun_merge)

## [1] "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"

proj4string(pro.sp)

## [1] "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"

proj4string(mun_merge)==proj4string(pro.sp)

## [1] TRUE

3 Creating a quantile map for the variable ANC with the borders of the provinces

tm_shape(mun_merge)+
  tm_fill(col="ANC",style = "quantile", n=3,  title = "ANC")+
  tm_layout(frame = FALSE,legend.outside = TRUE)+
 tm_borders(lwd=0.01)+
   tm_shape(pro.sp)+
  tm_borders(lwd=3)

tmap_save(filename = "ANC.png", height=5, width=5)

4 Spatial Analysis

4.1 creating a spatial weight matrix

nb<- poly2nb(mun_merge)
listw<- nb2listw(nb, style = "W")

4.2 Computing the Global Moran’s I

globalMoran <- moran.test(mun_merge$ANC, listw, na.action =  na.exclude)
globalMoran

## 
##  Moran I test under randomisation
## 
## data:  mun_merge$ANC  
## weights: listw    
## 
## Moran I statistic standard deviate = 10.211, p-value <
## 0.00000000000000022
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.600578644      -0.008695652       0.003560323

4.3 data in the shapefile

head(st_drop_geometry(st_as_sf( mun_merge)))

head(mun_merge@data)

4.4 Creating a table with the Local Moran’s I

Also Adding a quadrant variable to the local Moran’s I table

lmoran<- cbind(mun_merge@data, localmoran(mun_merge$ANC, listw, p.adjust.method="none", adjust.x=TRUE))
#lmoran

# centers the local Moran's around the mean
lmoran$Ii <- lmoran$Ii - mean(lmoran$Ii, na.rm = TRUE) 
lmoran$lag.ANC<-  lag.listw(listw,lmoran$ANC, NAOK = TRUE)

# centers the variable of interest around its mean

lmoran$ANCs <- lmoran$ANC - mean(lmoran$ANC, na.rm = TRUE) 
lmoran$lag.ANC <- lmoran$lag.ANC - mean(lmoran$lag.ANC, na.rm = TRUE) 

signif <- 0.05
#lmoran


lmoran <- lmoran%>% 
  mutate(quadrant= ifelse(ANCs>0 & lag.ANC > 0, 1, 0)) %>% 
  mutate(quadrant= ifelse(ANCs<0 & lag.ANC < 0, 2, quadrant)) %>% 
 mutate(quadrant= ifelse(ANCs<0 & lag.ANC > 0, 3, quadrant)) %>% 
 mutate(quadrant= ifelse(ANCs>0 & lag.ANC < 0, 4, quadrant)) %>%   mutate(quadrant= ifelse(lmoran$`Pr(z > 0)` > signif, 0, quadrant)) %>% 
  mutate(quadrant2= ifelse(ANC_st>0 & lagANC_st > 0, 1, 0)) %>% 
mutate(quadrant2= ifelse(ANC_st<0 & lagANC_st < 0, 2, quadrant2)) %>% 
  mutate(quadrant2= ifelse(ANC_st<0 & lagANC_st > 0, 3, quadrant2)) %>% 
 mutate(quadrant2= ifelse(ANC_st>0 & lagANC_st < 0, 4, quadrant2)) %>% 
 mutate(quadrant2= ifelse(lmoran$LISA_PANC > signif, 0, quadrant2))
  
mun_merge_new<- merge(mun_merge, lmoran, by.x="GID_3", by.y="GID_3")

4.5 comparing the geoda p-values and the R p-values

head(mun_merge_new@data %>% 
  select(LISA_PANC.y,`Pr(z > 0)` ))

4.6 LISA MAPS (P-VALUE < 0.05)

# R p value
 breaks = c(0, 1, 2, 3, 4, 5) 
LISA1<-tm_shape(mun_merge_new) + tm_fill(col = "quadrant", breaks = breaks, palette=  c("white","red","blue",rgb(0,0,1,alpha=0.4),rgb(1,0,0,alpha=0.4)), labels = c("Not significant", "High-High","Low-Low","Low-High","High-Low"), title="")+
  tm_legend(text.size = 1)  +
 # tm_scale_bar(position = c("LEFT", "BOTTOM"),text.size = 1.0)+
  # tm_compass(type = "8star",   position = c("RIGHT", "BOTTOM"),      show.labels = 2,   text.size = 0.5)+
    tm_borders(alpha=.5) +
   tm_layout( frame = FALSE,  title = "LISA with the R p-values ")

# Geoda p Value

breaks = c(0, 1, 2, 3, 4, 5) 
LISA2<- tm_shape(mun_merge_new) + tm_fill(col = "quadrant2", breaks = breaks, palette=  c("white","red","blue",rgb(0,0,1,alpha=0.4),rgb(1,0,0,alpha=0.4)), labels = c("Not significant", "High-High","Low-Low","Low-High","High-Low"), title="")+
  tm_legend(text.size = 1)  +
 # tm_scale_bar(position = c("LEFT", "BOTTOM"),text.size = 1.0)+
  # tm_compass(type = "8star",   position = c("RIGHT", "BOTTOM"),      show.labels = 2,   text.size = 0.5)+
    tm_borders(alpha=.5) +
   tm_layout( frame = FALSE, title = "LISA with the GEODA p-values")

tmap_arrange(LISA1, LISA2)

5 creating a nicer map with the p values of GEODA

 tm_shape(mun_merge_new) + tm_fill(col = "quadrant2", breaks = breaks, palette=  c("white","red","blue",rgb(0,0,1,alpha=0.4),rgb(1,0,0,alpha=0.4)), labels = c("Not significant", "High-High","Low-Low","Low-High","High-Low"), title="")+
  tm_legend(text.size = 1)  +
 # tm_scale_bar(position = c("LEFT", "BOTTOM"),text.size = 1.0)+
  # tm_compass(type = "8star",   position = c("RIGHT", "BOTTOM"),      show.labels = 2,   text.size = 0.5)+
    tm_borders(lwd =.05) +
    tm_shape(pro.sp)+
  tm_borders(lwd=3)+
   tm_layout( frame = FALSE, title = "LISA (GEODA p-values)")

tmap_save(filename = "ANC_lisa.png", height=5, width=5)