#### Due Date: February 12, 2014 Total Points: 40
## Use the Columbus, OH crime data from
## http://myweb.fsu.edu/jelsner/data/columbus.zip.
## 1. Read the data into R and create a choropleth map of the CRIME variable
## using the spplot method. 2. Compute Moran's I for the CRIME variable. 3.
## Create a Moran's scatter plot using the CRIME variable. Label the axes and
## include the regression line. 4. Check the value of Moran's I by
## determining the slope of the regression line.
## 1. Read the data into R and create a choropleth map of the CRIME variable
## using the spplot method. open file
require(maptools)
## Loading required package: maptools
## Warning: package 'maptools' was built under R version 3.0.2
## Loading required package: sp
## Warning: package 'sp' was built under R version 3.0.2
## Checking rgeos availability: TRUE
require(rgdal)
## Loading required package: rgdal
## Warning: package 'rgdal' was built under R version 3.0.2
## rgdal: version: 0.8-14, (SVN revision 496)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 1.10.1, released 2013/08/26
## Path to GDAL shared files: C:/Users/ADVQuant05/Documents/R/win-library/3.0/rgdal/gdal
## GDAL does not use iconv for recoding strings.
## Loaded PROJ.4 runtime: Rel. 4.8.0, 6 March 2012, [PJ_VERSION: 480]
## Path to PROJ.4 shared files: C:/Users/ADVQuant05/Documents/R/win-library/3.0/rgdal/proj
download.file("http://myweb.fsu.edu/jelsner/data/columbus.zip", "columbus.zip",
mode = "wb")
unzip("columbus.zip")
CC = readOGR(dsn = ".", layer = "columbus")
## OGR data source with driver: ESRI Shapefile
## Source: ".", layer: "columbus"
## with 49 features and 20 fields
## Feature type: wkbPolygon with 2 dimensions
## check data
slotNames(CC)
## [1] "data" "polygons" "plotOrder" "bbox" "proj4string"
str(CC, max.level = 2)
## Formal class 'SpatialPolygonsDataFrame' [package "sp"] with 5 slots
## ..@ data :'data.frame': 49 obs. of 20 variables:
## ..@ polygons :List of 49
## ..@ plotOrder : int [1:49] 21 9 5 20 1 40 6 42 2 7 ...
## ..@ bbox : num [1:2, 1:2] 5.87 10.79 11.29 14.74
## .. ..- attr(*, "dimnames")=List of 2
## ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slots
head(CC@data)
## AREA PERIMETER COLUMBUS_ COLUMBUS_I POLYID NEIG HOVAL INC CRIME
## 0 0.30944 2.441 2 5 1 5 80.47 19.531 15.73
## 1 0.25933 2.237 3 1 2 1 44.57 21.232 18.80
## 2 0.19247 2.188 4 6 3 6 26.35 15.956 30.63
## 3 0.08384 1.428 5 2 4 2 33.20 4.477 32.39
## 4 0.48889 2.997 6 7 5 7 23.23 11.252 50.73
## 5 0.28308 2.336 7 8 6 8 28.75 16.029 26.07
## OPEN PLUMB DISCBD X Y NSA NSB EW CP THOUS NEIGNO
## 0 2.8507 0.2172 5.03 38.80 44.07 1 1 1 0 1000 1005
## 1 5.2967 0.3206 4.27 35.62 42.38 1 1 0 0 1000 1001
## 2 4.5346 0.3744 3.89 39.82 41.18 1 1 1 0 1000 1006
## 3 0.3944 1.1869 3.70 36.50 40.52 1 1 0 0 1000 1002
## 4 0.4057 0.6246 2.83 40.01 38.00 1 1 1 0 1000 1007
## 5 0.5631 0.2541 3.78 43.75 39.28 1 1 1 0 1000 1008
df = CC@data
## SSPLOT
plot(CC)
spplot(CC, "CRIME", pretty = TRUE)
require(RColorBrewer)
## Loading required package: RColorBrewer
range(CC$CRIME)
## [1] 0.1783 68.8920
rng = seq(0, 70, 10)
cls = brewer.pal(7, "Greens")
spplot(CC, "CRIME", col.regions = cls, at = rng, sub = "Crime Data in Oho")
## 2.Compute Moran's I for the CRIME variable.
## BUILD NB
require(spdep)
## Loading required package: spdep
## Warning: package 'spdep' was built under R version 3.0.2
## Loading required package: Matrix
## Loading required package: lattice
CC.nb = poly2nb(CC)
CC.wts = nb2listw(CC.nb)
## CALCULATE MORAN'S I
m = length(CC$CRIME)
s = Szero(CC.wts)
moran(CC$CRIME, CC.wts, n = m, S0 = s)
## $I
## [1] 0.5002
##
## $K
## [1] 2.226
## TEST FOR WEIGHTING STYLE
moran.test(CC$CRIME, nb2listw(CC.nb, style = "W"))$estimate[1]
## Moran I statistic
## 0.5002
moran.test(CC$CRIME, nb2listw(CC.nb, style = "B"))$estimate[1] # binary
## Moran I statistic
## 0.5155
moran.test(CC$CRIME, nb2listw(CC.nb, style = "S"))$estimate[1] # variance-stabilizing
## Moran I statistic
## 0.5057
## 3. Create a Moran's scatter plot using the CRIME variable. Label the axes
## and include the regression line.
require(ggplot2)
## Loading required package: ggplot2
## Warning: package 'ggplot2' was built under R version 3.0.2
Wsids = lag.listw(CC.wts, CC$CRIME)
dat = data.frame(CC = CC$CRIME, Wsids = Wsids)
ggplot(dat, aes(x = CC, y = Wsids)) + geom_point() + geom_smooth(method = "lm") +
xlab("Columbus Crime") + ylab("Spatial Lag of crime data")
## 4. Check the value of Moran's I by determining the slope of the regression
## line.
lm(Wsids ~ CC, data = dat)
##
## Call:
## lm(formula = Wsids ~ CC, data = dat)
##
## Coefficients:
## (Intercept) CC
## 17.5 0.5