1 Objective

  1. Run various measures of spatial correlation
  2. Evaluate both global and local measures of spatial autocorrelation
  3. Identify spatial clusters

2 Key focus

  1. Moran Statistics
  2. Getis-Ord approach (Hot-spot)

3 Preliminary material

  1. Moran’s I
  2. Interpretations of Moran’s I
  3. Moran Scatter Plot
  4. Geary’s C
  5. Local Spatial Autocorrelation Clusters
  6. LISA Principle
  7. Local Moran
  8. Local G Statistics
  9. Issues and Interpretation
library(tidyverse)
library(sf)
library(sp)
library(spdep)
library(rgdal) # Bindings for the Geospatial Data Abstraction Library
library(rgeos) # Interface to Geometry Engine - Open Source
library(tmap)
library(tmaptools)
library(spgwr)
library(grid)
library(gridExtra)
library(leaflet)

options(prompt="R> ", digits=4, scipen=999)

4 Import datasets

4.1 Non-spatial data

Census.Data <- read.csv("practicaldata.csv")
glimpse(Census.Data)
## Rows: 749
## Columns: 5
## $ OA            <chr> "E00004120", "E00004121", "E00004122", "E00004123", "E00…
## $ White_British <dbl> 42.36, 47.20, 40.68, 49.66, 51.14, 41.42, 48.54, 48.68, …
## $ Low_Occupancy <dbl> 6.2937, 5.9322, 2.9126, 0.9259, 2.0000, 3.9326, 5.5556, …
## $ Unemployed    <dbl> 1.8939, 2.6882, 1.2121, 2.8037, 3.8168, 3.8462, 4.5455, …
## $ Qualification <dbl> 73.63, 69.90, 67.58, 60.78, 65.99, 74.21, 62.45, 60.35, …

4.2 Spatial data

Output.Areas <- readOGR(".", "Camden_oa11")
## Warning in OGRSpatialRef(dsn, layer, morphFromESRI = morphFromESRI, dumpSRS
## = dumpSRS, : Discarded datum Ordnance_Survey_of_Great_Britain_1936 in Proj4
## definition: +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000
## +y_0=-100000 +ellps=airy +units=m +no_defs
## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/admin/Desktop/Linh Data Studio/Spatial Data/Spatial Modelling", layer: "Camden_oa11"
## with 749 features
## It has 1 fields
glimpse(Output.Areas)
## Formal class 'SpatialPolygonsDataFrame' [package "sp"] with 5 slots
##   ..@ data       :'data.frame':  749 obs. of  1 variable:
##   .. ..$ OA11CD: chr [1:749] "E00004527" "E00004525" "E00004522" "E00004287" ...
##   ..@ polygons   :List of 749
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. .. [list output truncated]
##   ..@ plotOrder  : int [1:749] 154 604 597 195 670 418 303 243 97 676 ...
##   ..@ bbox       : num [1:2, 1:2] 523954 180960 531555 187604
##   .. ..- attr(*, "dimnames")=List of 2
##   ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
##   ..$ comment: chr "TRUE"
Output.Areas2 <- read_sf("Camden_oa11.shp")
glimpse(Output.Areas2)
## Rows: 749
## Columns: 2
## $ OA11CD   <chr> "E00004527", "E00004525", "E00004522", "E00004287", "E0000420…
## $ geometry <POLYGON [m]> POLYGON ((530648 181230, 53..., POLYGON ((530511 1815…

5 Merge datasets

OA.Census <- merge(Output.Areas, Census.Data, by.x="OA11CD", by.y="OA")

transform into a sf object

OA.Census_sf <- st_as_sf(OA.Census)

head(OA.Census_sf, 5)
## Simple feature collection with 5 features and 5 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 524300 ymin: 181200 xmax: 530700 ymax: 185100
## Projected CRS: Transverse_Mercator
##        OA11CD White_British Low_Occupancy Unemployed Qualification
## 397 E00004527         48.29        12.745      7.512         35.81
## 395 E00004525         40.94        16.807      5.991         42.41
## 392 E00004522         44.16         8.547      2.116         56.48
## 160 E00004287         31.87        12.613      3.286         67.86
## 80  E00004206         56.45        19.685      7.983         31.75
##                           geometry
## 397 POLYGON ((530648 181230, 53...
## 395 POLYGON ((530511 181531, 53...
## 392 POLYGON ((530207 181434, 53...
## 160 POLYGON ((524355 185054, 52...
## 80  POLYGON ((528718 184565, 52...

6 Spatial distribution

tm_shape(OA.Census_sf) + 
  tm_fill("Qualification",
          palette = "Reds", 
          style = "quantile", 
          title = "% with a Qualification") +
  tm_borders(alpha=.4)

7 Neighbour structure

7.1 Find queen neighbors

neighbours <- poly2nb(OA.Census)
neighbours
## Neighbour list object:
## Number of regions: 749 
## Number of nonzero links: 4342 
## Percentage nonzero weights: 0.774 
## Average number of links: 5.797
neighbours_sf <- poly2nb(OA.Census_sf)
neighbours_sf
## Neighbour list object:
## Number of regions: 749 
## Number of nonzero links: 4342 
## Percentage nonzero weights: 0.774 
## Average number of links: 5.797

7.3 Find rook neighbors

neighbours2 <- poly2nb(OA.Census, queen = FALSE)
neighbours2
## Neighbour list object:
## Number of regions: 749 
## Number of nonzero links: 4176 
## Percentage nonzero weights: 0.7444 
## Average number of links: 5.575

7.4 Plot rook neighbors

plot(OA.Census, border = 'lightgrey')
plot(neighbours2, coordinates(OA.Census), add=TRUE, col='blue')

7.5 Plot queen vs rook

plot(OA.Census, border = 'lightgrey')
plot(neighbours, coordinates(OA.Census), add=TRUE, col='red')
plot(neighbours2, coordinates(OA.Census), add=TRUE, col='blue')

8 Global spatial autocorrelation

8.1 Define weights matrix

listw <- nb2listw(neighbours2) # sf file
listw
## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 749 
## Number of nonzero links: 4176 
## Percentage nonzero weights: 0.7444 
## Average number of links: 5.575 
## 
## Weights style: W 
## Weights constants summary:
##     n     nn  S0    S1   S2
## W 749 561001 749 285.4 3114

8.2 Global Moran test

The correlation score is between -1 and 1. Much like a correlation coefficient, 1 determines perfect positive spatial autocorrelation (so our data is clustered), 0 identifies the data is randomly distributed and -1 represents negative spatial autocorrelation (so dissimilar values are next to each other).

globalMoran <- moran.test(OA.Census$Qualification, listw)
globalMoran
## 
##  Moran I test under randomisation
## 
## data:  OA.Census$Qualification  
## weights: listw    
## 
## Moran I statistic standard deviate = 24, p-value <0.0000000000000002
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##         0.5448699        -0.0013369         0.0005056
globalMoran[["estimate"]][["Moran I statistic"]]
## [1] 0.5449
globalMoran['p.value']
## $p.value
## [1] 0.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000001188

The Moran I statistic is 0.54, we can, therefore, determine that there our qualification variable is positively autocorrelated in Camden. In other words, the data does spatially cluster. We can also consider the p-value as a measure of the statistical significance of the model.

9 Local spatial autocorrelation

9.1 Moran scatterplot

moran <- moran.plot(OA.Census$Qualification, listw = nb2listw(neighbours2, style = "W"))

9.2 Compute local Moran

local <- localmoran(x = OA.Census$Qualification, listw = nb2listw(neighbours2, style = "W"))

Localmoran statistics: Ii: local moran statistic E.Ii: expectation of local moran statistic Var.Ii: variance of local moran statistic Z.Ii: standard deviate of local moran statistic Pr(): p-value of local moran statistic

9.3 Plot local Moran

9.3.1 A map the local moran statistic (Ii)

A positive value for Ii indicates that the unit is surrounded by units with similar values.

# binds results to our polygon shapefile
moran.map <- cbind(OA.Census, local)
tm_shape(moran.map) +
  tm_fill(col = "Ii",
         style = "quantile",
         title = "local moran statistic")

A positive value for Ii indicates that the unit is surrounded by units with similar values.

9.3.2 Plot LISA clusters

quadrant <- vector(mode = "numeric", length = nrow(local))

# centers the variables of interest around its mean
m.qualification <- OA.Census$Qualification - mean(OA.Census$Qualification)

# centres the local Moran's around the mean
m.local <- local[,1] - mean(local[,1])

# significant threshold
signif <- 0.1

# builds a data quadrant
quadrant[m.qualification >0 & m.local>0] <- 4  
quadrant[m.qualification <0 & m.local<0] <- 1      
quadrant[m.qualification <0 & m.local>0] <- 2
quadrant[m.qualification >0 & m.local<0] <- 3
quadrant[local[,5]>signif] <- 0   

# plot in r
brks <- c(0,1,2,3,4)
colors <- c("white","blue",rgb(0,0,1,alpha=0.4),rgb(1,0,0,alpha=0.4),"red")
plot(OA.Census,border="lightgray",col=colors[findInterval(quadrant,brks,all.inside=FALSE)])
box()
legend("bottomleft", legend = c("insignificant","low-low","low-high","high-low","high-high"),
       fill=colors,bty="n")

10 Getis-Ord Approach (Hot-spot)

The Getis-Ord Gi Statistic looks at neighbours within a defined proximity to identify where either high or low values cluster spatially. Here statistically significant hot-spots are recognised as areas of high values where other areas within a neighbourhood range also share high values too.

10.1 Search Radius

nb <- dnearneigh(coordinates(OA.Census), 0, 800) # 800 is radius

10.2 Create listw

nb_lw <- nb2listw(nb, style = "B")

10.3 Plot data and neighbours

plot(OA.Census, border = 'lightgrey')
plot(nb, coordinates(OA.Census), add=TRUE, col = 'red')

10.4 Getis-Ord Gi Statistic

local_g <- localG(OA.Census$Qualification, nb_lw)
local_g <- cbind(OA.Census, as.matrix(local_g))
names(local_g)[6] <- "gstat"

head(local_g, 5)
## class       : SpatialPolygonsDataFrame 
## features    : 5 
## extent      : 524326, 530660, 181181, 185111  (xmin, xmax, ymin, ymax)
## crs         : +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +units=m +no_defs 
## variables   : 6
## names       :    OA11CD,    White_British,    Low_Occupancy,       Unemployed,    Qualification,             gstat 
## min values  : E00004206, 31.8681318681319, 8.54700854700855, 2.11640211640212, 31.7460317460317, -3.12149898981916 
## max values  : E00004527, 56.4516129032258, 19.6850393700787, 7.98319327731092, 67.8571428571429,  1.18814041368312

10.5 Cluster map

tm_shape(local_g) + 
  tm_fill("gstat", 
          palette = "RdBu",
          style = "pretty") +
  tm_borders(alpha=.4)

The Gi Statistic is represented as a Z-score. Greater values represent a greater intensity of clustering and the direction (positive or negative) indicates high or low clusters. The final map should indicate the location of hot-spots across Camden. Repeat this for another variable.