date()
## [1] "Wed Feb 13 14:05:31 2013"
Use the Columbus, OH crime data from “http://geodacenter.org/downloads/data-files/columbus.zip”.
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")
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).
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")
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.