Exploratory Spatial Data Analysis (ESDA) with R

Juan Tomás Sayago Gómez

Análisis espacial’
Universidad Icesi - Cali
April 23, 2018

Juan Sayago ESDA Class 2018-01 1 / 82

1 Introduction to ESDA

2 Spatial Autocorrelation
Issues to consider
Hypothesis
Tests for the presence of Spatial autocorrelation
Spatial weights matrices
Global Spatial Autocorrelation
Moran’s I
Geary’s C
Moran’s Scatterplot
Local Spatial Autocorrelation

3 ESDA using R

Juan Sayago ESDA Class 2018-01 2 / 82

 Introduction to ESDA

What is Exploratory Spatial Data Analysis?

“ESDA is a collection of techniques to describe and visualize spatial

distributions, identify atypical locations or spatial outliers, discover patterns
of spatial association, clusters or hot-spots, and suggest spatial regimes or
other forms of spatial heterogeneity. Central to this conceptualization is
the notion of spatial autocorrelation or spatial association, i.e., the
phenomenon where locational similarity (observations in spatial proximity)
is matched by value similarity (attribute correlation)” Anselin(1998,79-80)

Juan Sayago ESDA Class 2018-01 3 / 82

 Introduction to ESDA

Introduction to Exploratory Data Analysis (EDA)

Graphical and visual methods or tools that are used to identify data
properties for purposes of
Detect patterns in the data
Formulate hypothesis from the data
Assess various aspects of models (e.g., goodness-of-fit)

Juan Sayago ESDA Class 2018-01 4 / 82

 Introduction to ESDA

Introduction to Exploratory Spatial Data Analysis (ESDA)

In many instances it is important to be able to link numerical and

graphical procedures with a map to answer questions such as ”Where are
those cases?”

Juan Sayago ESDA Class 2018-01 5 / 82

 Introduction to ESDA

Introduction to Exploratory Data Analysis (EDA)

EDA emphasizes on the interaction between human cognition and

computation in the form of dynamic statistical graphics that allow the
user to manipulate “views” of the data (via box-plots, scatterplots, etc).

EDA also uses statistics to identify outliers and detect patterns and other
characteristics in the data.

Juan Sayago ESDA Class 2018-01 6 / 82

 Introduction to ESDA

Introduction to Exploratory Spatial Data Analysis (ESDA)

It is an extension of EDA to detect spatial properties of data.

There is a need for additional techniques to

detect spatial patterns in data
formulate hypotheses based on the geography of the data
assess spatial models

Juan Sayago ESDA Class 2018-01 7 / 82

 Introduction to ESDA

Introduction to Exploratory Spatial Data Analysis (ESDA)

First example of ESDA by Mark

Monmonier (1989) Interactive
Spatial Data Analysis

With modern graphical interfaces

this is often done by “brushing” or
“highlighting” observations

Juan Sayago ESDA Class 2018-01 8 / 82

 Introduction to ESDA

Introduction to Exploratory Spatial Data Analysis (ESDA)

These geographic brushing tools are

now found in some GIS packages.
In addition to choropleth maps,
other displays/views include
histograms, box-plots, scatterplots
(see GeoDa or the Micromaps
package).

Juan Sayago ESDA Class 2018-01 9 / 82

 Introduction to ESDA

Introduction to Exploratory Spatial Data Analysis (ESDA)

“True ESDA pays attention to both spatial and attribute association”

(Anselin, 1998, p.79).

Juan Sayago ESDA Class 2018-01 10 / 82

 Spatial Autocorrelation

Spatial Autocorrelation

Tobler’s First Law of Geography

“All places are related but nearby places are more related than distant
place”
Spatial autocorrelation is the formal property that measures the degree to
which near and distant things are related

Juan Sayago ESDA Class 2018-01 11 / 82

 Spatial Autocorrelation Issues to consider

Issues to consider

We remember the types of data

Point data:
Accuracy of location is very important
Area/lattice dat:
Data reported for some regular or irregular areal unit

Juan Sayago ESDA Class 2018-01 12 / 82

 Spatial Autocorrelation Issues to consider

Issues to consider

There are 2 key components of spatial data:

Attribute data
Spatial data
And How you choose to sample, or aggregate, your data is very important

Affect spatial relationships in your data

Juan Sayago ESDA Class 2018-01 13 / 82

 Spatial Autocorrelation Issues to consider

Issues to consider

There are problems/challenges with areal data

Modifiable areal unit problem (MAUP)
Scale effect spatial data analysis at different scales may produce
different results
Zoning effect regrouping zones at a given scale may produce different
results
What are the ?
Optimal neighborhood/area size.
Alternative zoning schemes.

Juan Sayago ESDA Class 2018-01 14 / 82

 Spatial Autocorrelation Issues to consider

Issues to consider

Spatial autocorrelation is Scale-Dependent

Juan Sayago ESDA Class 2018-01 15 / 82

 Spatial Autocorrelation Hypothesis

Spatial Autocorrelation

Statistical test of match between locational similarity and attribute

similarity
Correlation between
Positive
Negative
Random (Absense)

Juan Sayago ESDA Class 2018-01 16 / 82

 Spatial Autocorrelation Hypothesis

Spatial Autocorrelation

Let’s remember from statistics class.

Juan Sayago ESDA Class 2018-01 17 / 82

 Spatial Autocorrelation Hypothesis

Spatial Autocorrelation

We need to consider how spatially correlated variables can be:

Which is negative? Which is positive? Which is random?

Juan Sayago ESDA Class 2018-01 18 / 82

 Spatial Autocorrelation Hypothesis

Spatial Autocorrelation

Let’s remember statistics again:

Null hypothesis:
Spatial randomness
Values observed at one location do not depend on values observed at
neighboring locations
Observed spatial pattern of values is equally likely as any other spatial
pattern
The location of values may be altered without affecting the
information content of the data

Juan Sayago ESDA Class 2018-01 19 / 82

 Spatial Autocorrelation Hypothesis

Spatial Autocorrelation

Positive spatial Autocorrelation

Clustering:
like values tend to be in similar locations
Neighbors are Similar
more alike than they would be under spatial randomness
Compatible with Diffusion:
but not necessarily caused by diffusion

Juan Sayago ESDA Class 2018-01 20 / 82

 Spatial Autocorrelation Hypothesis

Spatial Autocorrelation

Positive spatial Autocorrelation

Juan Sayago ESDA Class 2018-01 21 / 82

 Spatial Autocorrelation Hypothesis

Spatial Autocorrelation

Negative spatial Autocorrelation

Checkerboard Pattern
“opposite” of clustering
Neighbors are Dissimilar, but more dissimilar than they would be
under spatial randomness
Compatible with Competition, but not necessarily competition

Juan Sayago ESDA Class 2018-01 22 / 82

 Spatial Autocorrelation Hypothesis

Spatial Autocorrelation

Negative spatial Autocorrelation

Juan Sayago ESDA Class 2018-01 23 / 82

 Spatial Autocorrelation Hypothesis

Spatial Autocorrelation

Positive spatial Autocorrelation

Positive Spatial Autocorrelation Does Not Imply Diffusion
diffusion tends to yield positive spatial autocorrelation, but the
reverse is not necessary
spatial correlation may be due to structural factors, without contagion
or diffusion

Juan Sayago ESDA Class 2018-01 24 / 82

 Spatial Autocorrelation Hypothesis

Spatial Autocorrelation

Positive spatial Autocorrelation

What is the Cause Behind Clustering
True Contagion is the result of a contagious process, social
interaction, a dynamic process.
Apparent Contagion could be due to spatial heterogeneity,
stratification
It Cannot Be Distinguished In a Pure Cross Section

Juan Sayago ESDA Class 2018-01 25 / 82

 Spatial Autocorrelation Tests for the presence of Spatial autocorrelation

Global Spatial Autocorrelation Tests

Measure Degree of Clustering

Global Characteristic of Spatial Pattern - NOT Local
are like values more grouped in space than random
property of overall pattern = all the observations
test by means of a global spatial autocorrelation statistic
One statistic to summarize pattern in whole study area

Juan Sayago ESDA Class 2018-01 26 / 82

 Spatial Autocorrelation Tests for the presence of Spatial autocorrelation

Local Spatial Autocorrelation

Measure presence of Clusters or Hot-spots

Local Characteristic of Spatial Pattern
where are like values more grouped in space than random
property of local pattern = location-specific
test by means of a local spatial autocorrelation statistic
local clusters may be compatible with global spatial randomness
Location-specific statistics

Juan Sayago ESDA Class 2018-01 27 / 82

 Spatial Autocorrelation Spatial weights matrices

Spatial weights matrices

Neighborhoods can be defined in a number of ways Contiguity

(common boundary)
What is a “shared” boundary?
Distance (distance band, K-nearest neighbors)
How many “neighbors” to include, what distance do we use? General
weights (social distance, distance decay)
Other measures: Economic, social, political, trade.

Juan Sayago ESDA Class 2018-01 28 / 82

 Spatial Autocorrelation Spatial weights matrices

Spatial weights matrices

Juan Sayago ESDA Class 2018-01 29 / 82

 Spatial Autocorrelation Spatial weights matrices

Spatial weights matrices

Most common is using binary connectivity based on contiguity

wij = 1 if regions i and j are contiguous,
wij = 0 otherwise
May also be defined as a function of the distance between i and j
Distance of the line connecting the centroids of two areas
wij = dij−β
wij = exp (−βdij )

Juan Sayago ESDA Class 2018-01 30 / 82

 Spatial Autocorrelation Spatial weights matrices

Spatial weights matrices

Juan Sayago ESDA Class 2018-01 31 / 82

 Spatial Autocorrelation Spatial weights matrices

Spatial weights matrices

Juan Sayago ESDA Class 2018-01 32 / 82

 Spatial Autocorrelation Spatial weights matrices

Spatial weights matrices

A binary weights matrix looks like:
 
0 1 0 0
0 0 1 1
W = 1 1 0 0
 (1)
0 1 1 1
Observation 1 has neighbor 2
Observation 2 has neighbors 3 and 4
Observation 3 has neighbors 1 and 2
Observation 4 has neighbor 2, 3 and 4
A row-standardized matrix it looks like:
 
0 1 0 0
0 0 0.5 0.5
W = 0.5 0.5 0
 (2)
0
0 .33 .33 .33
Juan Sayago ESDA Class 2018-01 33 / 82
 Spatial Autocorrelation Global Spatial Autocorrelation

Spatial Autocorrelation Statistics

Formal Test of Match Between Value Similarity and Locational

Similarity
Statistic Summarizes Both Aspects
Significance
how likely is it (p-value) that the computed statistic would take this
(extreme) value in a spatially random pattern

Juan Sayago ESDA Class 2018-01 34 / 82

 Spatial Autocorrelation Global Spatial Autocorrelation

Spatial Autocorrelation Statistics

Tests for the presence of spatial autocorrelation

Global
Moran’s I
Geary’s C
Local (LISA – Local Indicators of Spatial Autocorrelation)
Local Moran’s
Getis Gi∗
Other tests that are more simple:
The Chi-square Test for Spatial Independence
The Join Count Statistic

Juan Sayago ESDA Class 2018-01 35 / 82

 Spatial Autocorrelation Global Spatial Autocorrelation

Moran’s Index

Classic/best measure of spatial autocorrelation

Depends upon definition of neighboring unit via the spatial weights
matrix
Typically ranges from −1 to 1
Like regression, it has a few assumptions
Regional x/y values all come from normal distributions w/same mean
and variance for each region
Randomly rearrange the data on map and compute I many times,
would have a normal distribution
Why? Because we use the normal distribution to calculate the p-value

Juan Sayago ESDA Class 2018-01 36 / 82

 Spatial Autocorrelation Global Spatial Autocorrelation

Moran’s Index

Pn Pn
n i=1 j=1 wij (yi − ȳ )(yj − ȳ )
I = Pn (3)
( i=1 (yi − ȳ )2 )(
PP
i6=j wij )

n = number of regions.
Wij = Measure of Spatial proximity between region i and j

Juan Sayago ESDA Class 2018-01 37 / 82

 Spatial Autocorrelation Global Spatial Autocorrelation

Moran’s Index

Juan Sayago ESDA Class 2018-01 38 / 82

 Spatial Autocorrelation Global Spatial Autocorrelation

Moran’s Index

Inference
Normal: assume uncorrelated normal distribution
Z-Value: compute E (I ) and Var (I )
z = (I − E (I ))/SD(I ) Regular Standarization.

Juan Sayago ESDA Class 2018-01 39 / 82

 Spatial Autocorrelation Global Spatial Autocorrelation

Moran’s Index

For Significant Statistics Only

Use z-value
The result Index depends on W
Positive S.A.:
zI > 0 for p < 0.05
no distinction between clustering of high or low values
Negative S.A.:
zI < 0 for p < 0.05

Juan Sayago ESDA Class 2018-01 40 / 82

 Spatial Autocorrelation Global Spatial Autocorrelation

Moran’s Index

Variance Instability of Rates

non-constant variance violates assumption of stationarity
may lead to spurious indication of spatial autocorrelation
Empirical Bayes Adjustment (Assuncao-Reis)
standardize each rate
use standardized rate

Juan Sayago ESDA Class 2018-01 41 / 82

 Spatial Autocorrelation Global Spatial Autocorrelation

Moran’s Index

Juan Sayago ESDA Class 2018-01 42 / 82

 Spatial Autocorrelation Global Spatial Autocorrelation

Geary’s Coefficient

Geary’s C typically ranges from 0 to 3

It Cannot be Negative
A Spatially uncorrelated process has an expected C = 1
Values less than 1 indicate positive spatial autocorrelation
Values greater than 1 indicate negative autocorrelation

Juan Sayago ESDA Class 2018-01 43 / 82

 Spatial Autocorrelation Global Spatial Autocorrelation

Geary’s Coefficient

Juan Sayago ESDA Class 2018-01 44 / 82

 Spatial Autocorrelation Global Spatial Autocorrelation

Geary’s Coefficient

Relationship of Moran’s I and Geary’s C

C approaches 0 and I approaches 1 when similar values are clustered
C approaches 3 and I approaches -1 when dissimilar values tend to
cluster
High values of C measures correspond to low values of I
So the two measures are inversely related

Juan Sayago ESDA Class 2018-01 45 / 82

 Spatial Autocorrelation Global Spatial Autocorrelation

Moran’s Scatterplot

Moran’s I as a Regression Slope

in matrix notation:
z 0 Wz
I = 0
zz
The Moran’s I is the slope in regression
Wz = a + I × z
Moran Scatter Plot
linear association between Wz on the y-axis and z on the x-axis
each point is pair (zi , Wzi ), slope is Moran’s I

Juan Sayago ESDA Class 2018-01 46 / 82

 Spatial Autocorrelation Global Spatial Autocorrelation

Moran’s Scatterplot

Juan Sayago ESDA Class 2018-01 47 / 82

 Spatial Autocorrelation Global Spatial Autocorrelation

Moran’s Scatterplot

Link with Local Spatial Autocorrelation

Four Categories of SA
Positive Spatial Autocorrelation
high-high and low-low: spatial clusters
Negative Spatial Autocorrelation
high-low and low-high: spatial outliers
Only Suggestive
no suggestion of significance
relative to mean

Juan Sayago ESDA Class 2018-01 48 / 82

 Spatial Autocorrelation Local Spatial Autocorrelation

Local Indicators of Spatial Autocorrelation – LISA

Anselin (1995)
Local Spatial Statistic
indicate significant spatial autocorrelation for each location
Local-Global Relation
sum of LISA proportional to a corresponding global indicator of
spatial autocorrelation

Juan Sayago ESDA Class 2018-01 49 / 82

 Spatial Autocorrelation Local Spatial Autocorrelation

LISA- Local Moran’s I

Used to determine if local autocorrelation exists around each region

Index
Returns an Index value for each region
Often used after global Moran’s Index to see if:
The study area is homogeneous (local statistics similar across regions)
There are local outliers that contribute to a significant global statistic

Juan Sayago ESDA Class 2018-01 50 / 82

 Spatial Autocorrelation Local Spatial Autocorrelation

LISA- Local Moran’s I

n(yi − ȳ ) X
Ii = P 2
wij (yi − ȳ ) (4)
j (yi − ȳ ) j

Juan Sayago ESDA Class 2018-01 51 / 82

 ESDA using R

Load required library and Data

First We need to install the packages spdep and maptools

> [Link](c('sp','spdep','maptools'))
The downloaded binary packages are in
/var/folders/qq/d9kxh5bj3xs7hs77g3fhrzp80000gn/T//Rtmp
After installing the libraries we need to load them up
> library(spdep)
> library(maptools)

Juan Sayago ESDA Class 2018-01 52 / 82

 ESDA using R

Load the Data

> shape <- readShapePoly("[Link]",IDvar="NAME" )
> summary(shape)
Object of class SpatialPolygonsDataFrame
Coordinates:
min max
x -124.7631 -66.94989
y 24.5231 49.38436
Is projected: NA
proj4string : [NA]
Data attributes:
SP_ID STATEFP STATENS AFFGEOID GEOID
Alabama : 1 01 : 1 00068085: 1 0400000US01: 1 01 : 1
Arizona : 1 04 : 1 00294478: 1 0400000US04: 1 04 : 1
Arkansas : 1 05 : 1 00448508: 1 0400000US05: 1 05 : 1
California : 1 06 : 1 00481813: 1 0400000US06: 1 06 : 1
Colorado : 1 08 : 1 00606926: 1 0400000US08: 1 08 : 1
Connecticut: 1 09 : 1 00662849: 1 0400000US09: 1 09 : 1
(Other) :43 (Other):43 (Other) :43 (Other) :43 (Other):43
STUSPS NAME LSAD ALAND AWATER
AL : 1 Alabama : 1 00:49 Min. :1.584e+08 Min. :1.863e+07
AR : 1 Arizona : 1 1st Qu.:9.279e+10 1st Qu.:1.536e+09
AZ : 1 Arkansas : 1 Median :1.389e+11 Median :3.543e+09
CA : 1 California : 1 Mean :1.562e+11 Mean :8.744e+09
CO : 1 Colorado : 1 3rd Qu.:2.062e+11 3rd Qu.:8.504e+09
CT : 1 Connecticut: 1 Max. :6.766e+11 Max. :1.040e+11
(Other):43 (Other) :43
State Violent_Cr Murder Forcible_r
Alabama : 1 Min. : 119.4 Min. : 0.900 Min. :12.90
Arizona : 1 1st Qu.: 299.3 1st Qu.: 3.000 1st Qu.:26.90
Arkansas : 1 Median : 385.1 Median : 4.900 Median :32.70
California : 1 Mean : 432.1 Mean : 5.486 Mean :33.37
Colorado Juan
: 1Sayago
3rd Qu.: 512.3 3rd Qu.: 6.600ESDA3rdClass
Qu.:38.90 2018-01 53 / 82
 ESDA using R

Create the Weight matrix- Queen

> nb2=poly2nb(shape, queen=TRUE)

> nb2

Neighbour list object:

Number of regions: 49
Number of nonzero links: 218
Percentage nonzero weights: 9.07955
Average number of links: 4.44898

> plot(shape, border="grey60")

> plot(nb2, coordinates(shape), add=TRUE, pcv=".", lwd=2)

Juan Sayago ESDA Class 2018-01 54 / 82

 ESDA using R

Create the Weight matrix- Queen

● ● ●
●
●
● ● ● ●●
● ●
● ●
● ● ●●
●
● ●
● ● ●
● ● ● ●●
●
● ● ●
● ● ●
● ● ●
● ● ●
●
● ● ●
● ●
●

Juan Sayago ESDA Class 2018-01 55 / 82

 ESDA using R

Create the Weight matrix- Rook

> nb2=poly2nb(shape, queen=FALSE)

> nb2

Neighbour list object:

Number of regions: 49
Number of nonzero links: 214
Percentage nonzero weights: 8.912953
Average number of links: 4.367347

> plot(shape, border="grey60")

> plot(nb2, coordinates(shape), add=TRUE, pcv=".", lwd=2)

Juan Sayago ESDA Class 2018-01 56 / 82

 ESDA using R

Create the Weight matrix- Rook

● ● ●
●
●
● ● ● ●●
● ●
● ●
● ● ●●
●
● ●
● ● ●
● ● ● ●●
●
● ● ●
● ● ●
● ● ●
● ● ●
●
● ● ●
● ●
●

Juan Sayago ESDA Class 2018-01 57 / 82

 ESDA using R

Create the Weight matrix- Rook

> sids <- readShapePoly("[Link]", ID="FIPSNO")
> sids_nbr<-poly2nb(sids,queen=F)
> plot(sids, border="grey60")
> plot(sids_nbr, coordinates(sids), add=TRUE, pcv=".", lwd=2)

● ● ● ● ● ●
●
●● ●
● ● ● ●● ●
● ● ● ● ● ● ●
●● ● ● ●● ●
● ● ●
● ● ● ● ● ● ● ● ●
● ● ● ●
● ● ● ● ● ●
● ● ● ● ● ● ● ●
● ● ●
● ● ● ● ●
● ● ● ● ● ● ● ● ● ● ● ●
● ● ●
● ● ● ● ● ● ● ●
● ● ●
● ● ●
● ●
●

Juan Sayago ESDA Class 2018-01 58 / 82

 ESDA using R

Create the Weight matrix- Rook

● ● ● ● ● ●● ● ● ● ● ● ● ● ●●
● ● ● ● ●
● ● ● ●
● ● ● ● ● ●
● ● ●
● ● ● ● ● ● ● ●
● ● ● ● ● ●
● ● ● ● ●
● ● ● ● ●
● ● ● ● ● ● ●
● ● ● ● ● ● ● ● ● ● ● ●
● ●
● ● ● ● ●
● ● ● ● ● ● ●
● ● ●
● ● ●
● ●
●

Juan Sayago ESDA Class 2018-01 59 / 82

 ESDA using R

Create the Weight matrix - Queen

> sids_nbq<-poly2nb(sids,queen=T)
> plot(sids, border="grey60")
> plot(sids_nbq, coordinates(sids), add=TRUE, pcv=".", lwd=2)

● ● ● ● ● ●● ● ● ● ● ● ● ● ●●
● ● ● ● ●
● ● ● ●
● ● ● ● ● ●
● ● ●
● ● ● ● ● ● ● ●
● ● ● ● ● ●
● ● ● ● ●
● ● ● ● ●
● ● ● ● ● ● ●
● ● ● ● ● ● ● ● ● ● ● ●
● ●
● ● ● ● ●
● ● ● ● ● ● ●
● ● ●
● ● ●
● ●
●

Juan Sayago ESDA Class 2018-01 60 / 82

 ESDA using R

Create the Weight matrix- Comparison

> sids_nbq<-poly2nb(sids,queen=T)
> sids_nbr<-poly2nb(sids,queen=F)
> plot(sids, border="grey60")
> plot(sids_nbq, coordinates(sids), add=TRUE, pcv=".",col="red", lwd=2)
> plot(sids_nbr, coordinates(sids), add=TRUE, pcv=".",col="blue", lwd=2)

● ● ● ● ● ●● ● ● ● ● ● ● ● ●●
● ● ● ● ●
● ● ● ●
● ● ● ● ● ●
● ● ●
● ● ● ● ● ● ● ●
● ● ●
● ● ● ● ● ● ● ●
● ●
● ● ● ● ●
● ● ● ● ● ● ● ● ●
● ●
●
● ● ● ● ● ● ●
● ● ● ● ●
● ● ● ● ● ● ●
● ● ●
● ● ●
● ●
●

Juan Sayago ESDA Class 2018-01 61 / 82

 ESDA using R

Create the Weight matrix - 4 Knearest Neighbors

> coords <- coordinates(shape)
> [Link] <- knearneigh(coords, k=4)
> plot(shape, border="grey60")
> plot(knn2nb([Link]), coordinates(shape), add=TRUE, pcv=".", lwd=2)

● ● ●
●
●
● ● ● ●●
● ●
● ●
● ● ●●
●
● ●
● ● ●
● ● ● ●●
●
● ● ●
● ● ●
● ● ●
● ● ●
●
● ● ●
● ●
●

Juan Sayago ESDA Class 2018-01 62 / 82

 ESDA using R

Create the Weight matrix - 3 Knearest Neighbors

> coords <- coordinates(shape)
> [Link] <- knearneigh(coords, k=3)
> plot(shape, border="grey60")
> plot(knn2nb([Link]), coordinates(shape), add=TRUE, pcv=".", lwd=2)

● ● ●
●
●
● ● ● ●●
● ●
● ●
● ● ●●
●
● ●
● ● ●
● ● ● ●●
●
● ● ●
● ● ●
● ● ●
● ● ●
●
● ● ●
● ●
●

Juan Sayago ESDA Class 2018-01 63 / 82

 ESDA using R

Binary W matrix or Weighted W matrix

We can standardize the weight matrix by row.

> nb2_B=nb2listw(nb2, style="B", [Link]=TRUE)
> nb2_W=nb2listw(nb2, style="W", [Link]=TRUE)
You can save the weights.
> [Link](nb2,"[Link]")

Juan Sayago ESDA Class 2018-01 64 / 82

 ESDA using R

Moran’s I and Moran’s Scatterplot

> [Link](shape$Murder, listw=nb2_W,alternative="[Link]")

Moran I test under randomisation

data: shape$Murder
weights: nb2_W

Moran I statistic standard deviate = 3.2206, p-value = 0.001279

alternative hypothesis: [Link]
sample estimates:
Moran I statistic Expectation Variance
0.207123212 -0.020833333 0.005009799

Juan Sayago ESDA Class 2018-01 65 / 82

 ESDA using R

Moran’s I and Moran’s Scatterplot

> moran=[Link](shape$Murder, listw=nb2_B)

> moran

Moran I test under randomisation

data: shape$Murder
weights: nb2_B

Moran I statistic standard deviate = 2.7617, p-value = 0.002875

alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.163703437 -0.020833333 0.004464824

Juan Sayago ESDA Class 2018-01 66 / 82

 ESDA using R

Moran’s I and Moran’s Scatterplot

Montecarlo Moran’s I
> [Link](1234)
> bperm=[Link](shape$Murder, listw=nb2_W, nsim=999)
> bperm

Monte-Carlo simulation of Moran I

data: shape$Murder
weights: nb2_W
number of simulations + 1: 1000

statistic = 0.20712, observed rank = 998, p-value = 0.002

alternative hypothesis: greater

Juan Sayago ESDA Class 2018-01 67 / 82

 ESDA using R

Geary’s Coefficient

> gearyR=[Link](shape$Murder, listw=nb2_W)

> gearyR

Geary C test under randomisation

data: shape$Murder
weights: nb2_W

Geary C statistic standard deviate = 3.1328, p-value = 0.0008658

alternative hypothesis: Expectation greater than statistic
sample estimates:
Geary C statistic Expectation Variance
0.55757044 1.00000000 0.01994487

Juan Sayago ESDA Class 2018-01 68 / 82

 ESDA using R

Moran’s I and Moran’s Scatterplot

> [Link](shape$Murder, listw=nb2_W,label=[Link](shape$STUSPS),xlab="Murder Rate", ylab="Spatially Lagge
>

VA ● MD ●
10
Spatially Lagged Murder Rate

●
●
8

●
●
●
●
● ● DC ●
●
●
●
● ●
6

●
● ●
●
● ● ● ●●
●
●
● ●
●
● ● ● ●
4

●
●
●
● ● ● ●
●
● ●● ●
●
2

0 5 10 15 20 25 30

Murder Rate

Juan Sayago ESDA Class 2018-01 69 / 82

 ESDA using R

Moran’s I and Moran’s Scatterplot

> moran.plot2 <- function(x,wfile)
+ {
+ xname <- deparse(substitute(x)) # get name of variable
+ zx <- (x - mean(x))/sd(x)
+ wzx <- [Link](wfile,zx, [Link]=TRUE)
+ morlm <- lm(wzx ~ zx)
+ aa <- morlm$coefficients[1]
+ mori <- morlm$coefficients[2]
+ par(pty="s")
+ plot(zx,wzx,xlab=xname,ylab=paste("Spatial Lag of ",xname))
+ abline(aa,mori,col=2)
+ abline(h=0,lty=2,col=4)
+ abline(v=0,lty=2,col=4)
+ title(paste("Moran Scatterplot I= ",format(round(mori,4))))
+ }
> moran.plot2(shape$Murder, nb2_W)

Moran Scatterplot I= 0.2071

● ●
1.0
Spatial Lag of shape$Murder

●
●
●
0.5

●
●
●● ● ●
●
●
●
●●
●
0.0

● ●
● ●
● ● ● ●●
●
● ●
●
●● ●●
●
●
−0.5

●
●●
Juan Sayago
●●● ESDA Class 2018-01 70 / 82
 ESDA using R

Moran’s I and Moran’s Scatterplot

> moran.plot2(shape$Murder, nb2_W)

Moran Scatterplot I= 0.2071

● ●
1.0
Spatial Lag of shape$Murder

●
●
●
0.5

●
●
●● ● ●
●
●
●
●●
●
0.0

● ●
● ●
● ● ● ●●
●
● ●
●
●● ●●
●
●
−0.5

●
●● ●●
●
● ●● ●
●
●
−1.0

−1 0 1 2 3 4 5 6

shape$Murder

Juan Sayago ESDA Class 2018-01 71 / 82

 ESDA using R

Moran’s I and Moran’s Scatterplot

> moran.plot3 <- function(x,y,wfile)
+ {
+ xname <- deparse(substitute(x)) # get name of variable
+ yname =deparse(substitute(y))
+ zx <- (x - mean(x))/sd(x)
+ zy =(y - mean(y))/sd(y)
+ wzy <- [Link](wfile,zy, [Link]=TRUE)
+ morlm <- lm(wzy ~ zx)
+ aa <- morlm$coefficients[1]
+ mori <- morlm$coefficients[2]
+ par(pty="s")
+ plot(zx,wzy,xlab=xname,ylab=paste("Spatial Lag of ",yname))
+ abline(aa,mori,col=2)
+ abline(h=0,lty=2,col=4)
+ abline(v=0,lty=2,col=4)
+ title(paste("Moran Scatterplot I= ",format(round(mori,4))))
+ }
> moran.plot3(shape$Murder,shape$Robbery, nb2_W)

Moran Scatterplot I= 0.1859

●
●
1.0
l Lag of shape$Robbery

0.5

●●
● ●●
● ●
● ●
●
● ●
● ●
●
●
● ●●
●
0.0

●
● ●
●
● ●●
●
●
● ●
Juan Sayago
●
● ●
● ESDA Class 2018-01 72 / 82
 ESDA using R

Moran’s I and Moran’s Scatterplot

> moran.plot3(shape$Murder,shape$Robbery, nb2_W)

Moran Scatterplot I= 0.1859

●
●
1.0
Spatial Lag of shape$Robbery

0.5

●●
● ●●
● ●
● ●
●
● ●
● ●
●
●
● ●●
●
0.0

●
● ●
●
● ●●
●
●
● ●
● ●
● ●
● ●
−0.5

●
●
●●
● ● ●
●
−1.0

−1 0 1 2 3 4 5 6

shape$Murder

Juan Sayago ESDA Class 2018-01 73 / 82

 ESDA using R

Local Indicators of Spatial Association

> [Link] <- function(var,listw,signif,mapa) {

+ [Link] <- localmoran(var,listw, [Link]=T)
+ [Link] <- var - mean(var)
+ [Link] <- [Link][,1] - mean([Link][,1])
+ quadrant <- vector(mode="numeric",length=nrow([Link]))
+ wzx <- [Link](listw,var, [Link]=TRUE)
+ c.mI1 <- wzx - mean(wzx)
+ quadrant[[Link]>0 & c.mI1>0] <- 4
+ quadrant[[Link]<0 & c.mI1<0] <- 1
+ quadrant[[Link]<0 & c.mI1>0] <- 2
+ quadrant[[Link]>0 & c.mI1<0] <- 3
+ quadrant[[Link][,5]>signif] <- 0
+ 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(mapa,border="gray90",col=colors[findInterval(quadrant,brks,[Link]=FALSE)])
+ box()
+ legend("bottomright",legend=c("not significant","low-low","low-high","high-low","high-high"), fill=colors,bt
+ title("LISA Mapa de Clusters")
+ }

Juan Sayago ESDA Class 2018-01 74 / 82

 ESDA using R

Local Indicators of Spatial Association

Using 0.1 significance threshold.

> [Link](shape$Murder, nb2_W,0.1, shape)

LISA Mapa de Clusters

not significant
low−low
low−high
high−low
high−high

Juan Sayago ESDA Class 2018-01 75 / 82

 ESDA using R

Local Indicators of Spatial Association

Using 0.05 significance threshold.

> [Link](shape$Murder, nb2_W,0.05, shape)

LISA Mapa de Clusters

not significant
low−low
low−high
high−low
high−high

Juan Sayago ESDA Class 2018-01 76 / 82

 ESDA using R

Load the Data - Example with counties

> shapec <- readShapePoly("[Link]",IDvar="GEOID" )
> summary(shapec)
Object of class SpatialPolygonsDataFrame
Coordinates:
min max
x -82.64459 -71.77749
y 37.20154 45.01586
Is projected: NA
proj4string : [NA]
Data attributes:
SP_ID STATEFP COUNTYFP COUNTYNS GEOID
34001 : 1 34:21 001 : 4 00882228: 1 34001 : 1
34003 : 1 36:62 003 : 4 00882229: 1 34003 : 1
34005 : 1 42:67 005 : 4 00882230: 1 34005 : 1
34007 : 1 54:55 007 : 4 00882231: 1 34007 : 1
34009 : 1 009 : 4 00882232: 1 34009 : 1
34011 : 1 011 : 4 00882233: 1 34011 : 1
(Other):199 (Other):181 (Other) :199 (Other):199
NAME NAMELSAD LSAD CLASSFP MTFCC
Jefferson: 3 Jefferson County: 3 06:205 H1:199 G4020:205
Mercer : 3 Mercer County : 3 H6: 6
Monroe : 3 Monroe County : 3
Warren : 3 Warren County : 3
Wayne : 3 Wayne County : 3
Wyoming : 3 Wyoming County : 3
(Other) :187 (Other) :187
CSAFP CBSAFP METDIVFP FUNCSTAT ALAND
408 :32 35620 : 25 35614 : 14 A:199 Min. :5.868e+07
428 :13 37980 : 9 35084 : 7 C: 6 1st Qu.:9.234e+08
430 :11 38300 : 7 15804 : 3 Median :1.343e+09
104 :10 40380 : 6 33874 : 3 Mean :1.557e+09
170 : 8Juan10580
Sayago : 5 20524 : 2 ESDA
3rd Class
Qu.:2.070e+09 2018-01 77 / 82
 ESDA using R

Load the Data - Example with counties

Juan Sayago ESDA Class 2018-01 78 / 82

 ESDA using R

Create the Weight matrix- Queen

> nb2c=poly2nb(shapec, queen=TRUE)
> nb2c
Neighbour list object:
Number of regions: 205
Number of nonzero links: 1080
Percentage nonzero weights: 2.569899
Average number of links: 5.268293
> plot(shapec, border="grey60")
> plot(nb2c, coordinates(shapec), add=TRUE, pcv=".", lwd=2)
> nb2_Wc=nb2listw(nb2c, style="W", [Link]=TRUE)

●
● ●

● ●
● ●
● ●
● ● ● ● ● ●
● ● ●
● ● ● ● ●
● ● ● ●● ● ●
● ● ● ● ● ●
●
● ● ● ● ●● ● ● ● ● ●
●
● ● ● ● ● ● ●
● ● ● ●
● ● ● ● ● ● ● ● ●
● ● ● ● ● ●
● ● ●● ● ● ●●●
● ● ● ● ● ● ●● ●
● ● ● ● ● ●● ● ●
● ●
● ●
●
● ● ● ● ●● ●ESDA ● ● ● ●●●●
Juan Sayago ● ● ● Class ● ● ●● 2018-01 79 / 82
 ESDA using R

Moran’s I and Moran’s Scatterplot - Counties

> moran.plot2(shapec$Bachelor, nb2_Wc)

Moran Scatterplot I= 0.5038

●
2.0

●
●
●
●
1.5

● ● ●
●
Spatial Lag of shapec$Bachelor

● ● ●
●
●
●● ●
●
●● ●
1.0

●
●
● ●
● ●
●
● ● ●
● ● ● ● ●
●
●● ● ●
● ● ● ●
●
0.5

●● ●● ● ●
● ● ● ● ●
● ●
●● ●
●
● ● ●● ●●● ●
●● ● ● ● ●
●
● ●● ● ●
● ●
−1.5 −1.0 −0.5 0.0

● ● ● ●
● ● ● ●
● ●● ● ● ●
●●● ● ●● ●
● ● ● ● ●● ●
● ●● ● ● ● ● ● ●
●
●●
● ● ● ● ●
●●●
●● ● ●●● ● ●
● ● ●●● ●
● ●
● ● ● ●● ●● ● ●
●●●
●● ●●●● ●●
● ● ●● ● ● ●●
● ● ●●
● ● ●● ● ● ●
● ● ●
● ● ●
●
●
●
● ●
●

−2 −1 0 1 2 3

shapec$Bachelor

Juan Sayago ESDA Class 2018-01 80 / 82

 ESDA using R

LISA map counties

Using 0.05 significance threshold.

> [Link](shapec$Bachelor, nb2_Wc,0.05, shapec)

LISA Mapa de Clusters

not significant
low−low
low−high
high−low
high−high

Juan Sayago ESDA Class 2018-01 81 / 82

 ESDA using R

LISA map counties

Taller hacer Moran y LISA de NBI y Cobertura de energia electrica en

Colombia
1 Hacer Indice de Moran, diagrama de dispersión de Moran y LISA para
el NBI de Colombia
2 Usando los datos de desplazamiento en el shapefile
“shapec [Link]” hacer un mapa de expulsion, estandarizar la
variable, y hacer análisis de autocorrelación espacial.
3 Elegir una variable y hacer repetir el análisis.

Juan Sayago ESDA Class 2018-01 82 / 82