Create a LISA map Exercise

load data and libraries

load("C:/Users/molly/Documents/1_Coursework/courses by topic/Geography/Methods in Geography/Data/soco.rda")
library(maptools)
## Warning: package 'maptools' was built under R version 2.15.3
## Loading required package: foreign
## Loading required package: sp
## Loading required package: grid
## Loading required package: lattice
## Checking rgeos availability: FALSE Note: when rgeos is not available,
## polygon geometry computations in maptools depend on gpclib, which has a
## restricted licence. It is disabled by default; to enable gpclib, type
## gpclibPermit()
library(rgdal)
## Warning: package 'rgdal' was built under R version 2.15.3
## rgdal: version: 0.8-5, (SVN revision 449) Geospatial Data Abstraction
## Library extensions to R successfully loaded Loaded GDAL runtime: GDAL
## 1.9.2, released 2012/10/08 Path to GDAL shared files: C:/Program
## Files/R/R-2.15.2/library/rgdal/gdal GDAL does not use iconv for recoding
## strings. Loaded PROJ.4 runtime: Rel. 4.7.1, 23 September 2009,
## [PJ_VERSION: 470] Path to PROJ.4 shared files: C:/Program
## Files/R/R-2.15.2/library/rgdal/proj
library(spdep)
## Warning: package 'spdep' was built under R version 2.15.3
## Loading required package: boot
## Attaching package: 'boot'
## The following object(s) are masked from 'package:lattice':
## 
## melanoma
## Loading required package: Matrix
## Loading required package: MASS
## Loading required package: nlme
## Loading required package: deldir
## deldir 0.0-21
## Loading required package: coda
## Warning: package 'coda' was built under R version 2.15.3
## Loading required package: splines
library(classInt)
## Warning: package 'classInt' was built under R version 2.15.3
## Loading required package: class
## Loading required package: e1071
## Warning: package 'e1071' was built under R version 2.15.3
library(spdep)
library(gstat)
## Warning: package 'gstat' was built under R version 2.15.3

prepare the data

# create row-standardized Queens contiguity weights matrix
soco_nbq <- poly2nb(soco)  #queen's neighborhood
soco_nbq_w <- nb2listw(soco_nbq)
locm <- localmoran(soco$PPOV, soco_nbq_w)  #calculate the local moran's I
summary(locm)
##        Ii              E.Ii               Var.Ii            Z.Ii       
##  Min.   :-1.778   Min.   :-0.000722   Min.   :0.0901   Min.   :-3.983  
##  1st Qu.: 0.012   1st Qu.:-0.000722   1st Qu.:0.1420   1st Qu.: 0.029  
##  Median : 0.219   Median :-0.000722   Median :0.1658   Median : 0.517  
##  Mean   : 0.589   Mean   :-0.000722   Mean   :0.1861   Mean   : 1.408  
##  3rd Qu.: 0.746   3rd Qu.:-0.000722   3rd Qu.:0.1990   3rd Qu.: 1.761  
##  Max.   : 9.335   Max.   :-0.000722   Max.   :0.9981   Max.   :18.977  
##    Pr(z > 0)     
##  Min.   :0.0000  
##  1st Qu.:0.0391  
##  Median :0.3025  
##  Mean   :0.2900  
##  3rd Qu.:0.4886  
##  Max.   :1.0000

# manually make a moran plot standarize variables
soco$sPPOV <- scale(soco$PPOV)  #save to a new column

# create a lagged variable
soco$lag_sPPOV <- lag.listw(soco_nbq_w, soco$sPPOV)

summary(soco$sPPOV)
##        V1        
##  Min.   :-2.151  
##  1st Qu.:-0.700  
##  Median :-0.123  
##  Mean   : 0.000  
##  3rd Qu.: 0.591  
##  Max.   : 4.062
summary(soco$lag_sPPOV)
##        V1         
##  Min.   :-1.9925  
##  1st Qu.:-0.5298  
##  Median :-0.0762  
##  Mean   : 0.0085  
##  3rd Qu.: 0.4514  
##  Max.   : 2.6478
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")

# check out the outliers click on one or twon and then hit escape (or
# click finish)
identify(soco$sPPOV, soco$lab_sPPOV, soco$CNTY_ST, cex = 0.8)

plot of chunk unnamed-chunk-2

## integer(0)

# identify the moran plot quadrant for each observation
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  #WE ASSIGN A 5 TO ALL NON-SIGNIFICANT OBSERVATIONS

MAP THE RESULTS (courtesy of Paul Voss)

# Set the breaks for the thematic map classes
breaks <- seq(1, 5, 1)

# Set the corresponding labels for the thematic map classes
labels <- c("high-High", "low-Low", "High-Low", "Low-High", "Not Signif.")

# see ?findInterval - This is necessary for making a map
np <- findInterval(soco$quad_sig, breaks)

# Assign colors to each map class
colors <- c("red", "blue", "lightpink", "skyblue2", "white")
plot(soco, col = colors[np])  #colors[np] manually sets the color for each county
mtext("Local Moran's I", cex = 1.5, side = 3, line = 1)
legend("topleft", legend = labels, fill = colors, bty = "n")

plot of chunk unnamed-chunk-3