Applied Spatial Statistics: Problem Set # 1

Loury Migliorelli

date()

## [1] "Wed Feb 13 09:15:48 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))
tmp = download.file("http://geodacenter.org/downloads/data-files/columbus.zip", 
    "columbus.zip", mode = "wb")
unzip("columbus.zip")
CC = readShapeSpatial("columbus/columbus", proj4string = CRS("+proj=longlat +ellps=clrk66"))
slotNames(CC)

## [1] "data"        "polygons"    "plotOrder"   "bbox"        "proj4string"

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

suppressMessages(install.packages("RColorBrewer"))

## Error: trying to use CRAN without setting a mirror

suppressMessages(require(RColorBrewer))
range(CC$CRIME)

## [1]  0.1783 68.8920

rng = seq(0, 70, 10)
cls = brewer.pal(7, "Blues")
spplot(CC, "CRIME", col.regions = cls, at = rng)

1. Compute Moran's I for the CRIME variable.

suppressMessages(require(spatstat))
suppressMessages(require(spdep))
CC.nb = poly2nb(CC)
CC.wts = nb2listw(CC.nb)
m = length(CC$CRIME)
s = Szero(CC.wts)
moran.test(CC$CRIME, CC.wts)

## 
##  Moran's I test under randomisation
## 
## data:  CC$CRIME  
## weights: CC.wts  
##  
## Moran I statistic standard deviate = 5.589, p-value = 1.139e-08
## alternative hypothesis: greater 
## sample estimates:
## Moran I statistic       Expectation          Variance 
##          0.500189         -0.020833          0.008689

The Moran's I for the CRIME variable is 0.5002.

1. Create a Moran's scatter plot using the CRIME variable. Label the axes and include the regression line.

suppressMessages(require(ggplot2))
Wcc = lag.listw(CC.wts, CC$CRIME)
Wcc[1]

## [1] 24.71

j = CC.wts$neighbours[[1]]
j

## [1] 2 3

sum(CC$CRIME[j])/length(j)

## [1] 24.71

dat = data.frame(cc = CC$CRIME, Wcc = Wcc)
ggplot(dat, aes(x = cc, y = Wcc)) + geom_point() + geom_smooth(method = "lm") + 
    xlab("CRIME") + ylab("Spatial Lag of CRIME")

1. Check the value of Moran's I by finding the slope of the regression line.

lm(Wcc ~ cc, data = dat)

## 
## Call:
## lm(formula = Wcc ~ cc, data = dat)
## 
## Coefficients:
## (Intercept)           cc  
##        17.5          0.5

The slope of the regression line is 0.5002 which is equivalent to the Moran's I calculated for the CRIME variable.