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)
## 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")