Applied Spatial Statistics: Problem Set # 1

Holly Widen

date()
## [1] "Wed Feb 13 14:05:31 2013"

Due Date: February 13, 2013

Total Points: 40

Use the Columbus, OH crime data from “http://geodacenter.org/downloads/data-files/columbus.zip”.

  1. Read the data into R and create a choropleth map of the CRIME variable using the ssplot method.
suppressMessages(require(maptools))
## Warning: package 'sp' was built under R version 2.15.2
tmp = download.file("http://geodacenter.org/downloads/data-files/columbus.zip", 
    "columbus.zip", mode = "wb")
unzip("columbus.zip")
Ccrime = readShapeSpatial("columbus")
slotNames(Ccrime)
## [1] "data"        "polygons"    "plotOrder"   "bbox"        "proj4string"
str(Ccrime, max.level = 2)
## Formal class 'SpatialPolygonsDataFrame' [package "sp"] with 5 slots
##   ..@ data       :'data.frame':  49 obs. of  20 variables:
##   .. ..- attr(*, "data_types")= chr [1:20] "N" "N" "N" "N" ...
##   ..@ 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(Ccrime@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
suppressMessages(require(RColorBrewer))
range(Ccrime$CRIME)
## [1]  0.1783 68.8920
rng = seq(0, 70, 10)
cls = brewer.pal(7, "Purples")
spplot(Ccrime, "CRIME", col.regions = cls, at = rng, colorkey = list(space = "bottom"), 
    sub = "Residential Burglaries and Vehicle Thefts per 1000 Households")

plot of chunk choromap

  1. Compute Moran's I for the CRIME variable.
suppressMessages(require(spdep))
## Warning: package 'spdep' was built under R version 2.15.2
## Warning: package 'deldir' was built under R version 2.15.2
## Warning: package 'coda' was built under R version 2.15.2
Ccrime.nb = poly2nb(Ccrime)
Ccrime.wts = nb2listw(Ccrime.nb)
m = length(Ccrime$CRIME)
s = Szero(Ccrime.wts)
moran(Ccrime$CRIME, Ccrime.wts, n = m, S0 = s)
## $I
## [1] 0.5002
## 
## $K
## [1] 2.226

Moran's I is about 0.50 (indicating fairly high spatial autocorrelation) and the sample kurtosis is approximately 2.23 (indicating that the normality assumption can be used to make inferences about I).

  1. Create a Moran's scatter plot using the CRIME variable. Label the axes and include the regression line.
crime = Ccrime$CRIME
SLcrime = lag.listw(Ccrime.wts, crime)
suppressMessages(require(ggplot2))
dat = data.frame(crime = crime, SLcrime = SLcrime)
ggplot(dat, aes(x = crime, y = SLcrime)) + geom_point() + geom_smooth(method = "lm") + 
    xlab("Crime per 1000 Households") + ylab("Spatial Lag of Crime per 1000 Households")

plot of chunk MoransScatter

  1. Check the value of Moran's I by finding the slope of the regression line.
lm(SLcrime ~ crime, data = dat)
## 
## Call:
## lm(formula = SLcrime ~ crime, data = dat)
## 
## Coefficients:
## (Intercept)        crime  
##        17.5          0.5

The slope of the regression line is 0.50 which is equal to the value of Moran's I.