Key Concepts covered in the lecture include:
1. Regression requires spatial dependence (no autocorrelation!)
2. Moran's I Significane Testing through simulation
3. Moran's Scatterplot
4. Local Indicators of Spatial Autocorrelation (LISA)
5. Correlogram
6. Semivariograms
7. The Spatial Lag Regression Model
The in-class exercise focuses on Local and Global Spatial Autocorrelation (Moran's I).
Part 1: Copute Global Moran's I and create a LISA map
library(spdep)
## Loading required package: sp
## Loading required package: boot
## Loading required package: Matrix
## Loading required package: lattice
## Attaching package: 'lattice'
## The following object(s) are masked from 'package:boot':
##
## melanoma
## Loading required package: MASS
## Loading required package: nlme
## Loading required package: maptools
## Loading required package: foreign
## Loading required package: grid
## Checking rgeos availability: TRUE
## Loading required package: deldir
## deldir 0.0-21
## Loading required package: coda
## Loading required package: splines
library(classInt)
## Loading required package: class
## Loading required package: e1071
library(gstat)
## Loading required package: spacetime
## Loading required package: xts
## Loading required package: zoo
## Attaching package: 'zoo'
## The following object(s) are masked from 'package:base':
##
## as.Date, as.Date.numeric
# Add data and create a weights matrix
load("E:/Quant/inClassExercises/InClassExerciseData/soco.rda")
soco_nbq <- poly2nb(soco) #Queens neighbors
soco_nbq_w <- nb2listw(soco_nbq) #row standardize
# This is a sumulated Moran's I
moran.mc(soco$PERPOV, listw = soco_nbq_w, nsim = 999)
##
## Monte-Carlo simulation of Moran's I
##
## data: soco$PERPOV
## weights: soco_nbq_w
## number of simulations + 1: 1000
##
## statistic = 0.4347, observed rank = 1000, p-value = 0.001
## alternative hypothesis: greater
# Assessing Moran's I
MyMoran <- moran.mc(soco$PPOV, listw = soco_nbq_w, nsim = 999)
# plot Moran's I
hist(MyMoran$res, breaks = 50)
MyMoran #myMoran is .589 - look how far it is from the simulated histogram!
##
## Monte-Carlo simulation of Moran's I
##
## data: soco$PPOV
## weights: soco_nbq_w
## number of simulations + 1: 1000
##
## statistic = 0.5893, observed rank = 1000, p-value = 0.001
## alternative hypothesis: greater
# Plotting Moran's I (looking for outliers...)
mp <- moran.plot(soco$PPOV, soco_nbq_w, labels = as.character(soco$CNTY_ST),
xlab = "Percent Poverty", ylab = "Lag of Percent Poverty")
# Can identify outliers (important observations) by creating a standard
# ellipse
# We can also create this manually to have a little more control over it.
soco$sPPOV <- scale(soco$PPOV)
soco$lag_sPPOV <- lag.listw(soco_nbq_w, soco$sPPOV)
plot(x = soco$sPPOV, y = soco$lag_sPPOV, main = "Moran Scatterplot PPOV")
abline(h = 0, v = 0)
abline(lm(soco$lag_sPPOV ~ soco$sPPOV), lty = 3, lwd = 4, col = "red")
# Can manually click on a point and label it. A nice interactive tool.
# identify(soco$sPPOV, soco$lag_sPPOV, soco$CNTY_ST, cex=0.8)
# Create a LISA map
locm <- localmoran(soco$PPOV, soco_nbq_w) #Calculate the local Morann
# Organize the data...
soco$quad_sig <- NA
soco@data[(soco$sPPOV >= 0 & soco$lag_sPPOV >= 0) & (locm[, 5] <= 0.05), "quad_sig"] <- 1
soco@data[(soco$sPPOV <= 0 & soco$lag_sPPOV <= 0) & (locm[, 5] <= 0.05), "quad_sig"] <- 2
soco@data[(soco$sPPOV >= 0 & soco$lag_sPPOV <= 0) & (locm[, 5] <= 0.05), "quad_sig"] <- 3
soco@data[(soco$sPPOV >= 0 & soco$lag_sPPOV <= 0) & (locm[, 5] <= 0.05), "quad_sig"] <- 4
soco@data[(soco$sPPOV <= 0 & soco$lag_sPPOV >= 0) & (locm[, 5] <= 0.05), "quad_sig"] <- 5
# Set up the map...
breaks <- seq(1, 5, 1)
labels <- c("high-High", "low-Low", "High-Low", "Low-High", "Not Signif.")
np <- findInterval(soco$quad_sig, breaks)
colors <- c("red", "blue", "lightpink", "skyblue2", "white")
plot(soco, col = colors[np])
mtext("Local Moran's I", cex = 1.5, side = 3, line = 1)
legend("topleft", legend = labels, fill = colors, bty = "n")
# NOTE: There are no H-L or L-H features plotted because they are not
# significant (though they seemed to be pretty far from the Moran's
# scatterplot)
Part 2: Create a variogram from the PPOV variable
# Note that these can be computational intense and will take a while to
# create
# Create a correlogram of Moran's I
cor10 <- sp.correlogram(soco_nbq, soco$PPOV, order = 10, method = "I")
# cor10
plot(cor10)
# One step further... a variogram
plot(variogram(soco$PPOV ~ 1, locations = coordinates(soco), data = soco, cloud = F),
type = "b")
# It has a sill at about 0.008; a nugget effect of about 0.004. The range
# is about 300,000.