Add p-values from Local Moran's test to a shapefile with counts.

localMorans.R

library("foreign")
library("rgdal")
library("sp")
library("spdep")
library("maptools")
library("SpatialEpi")
library("plyr")

data <- read.dbf("counts.dbf")
shapes <- readShapePoly("counts.shp")

#create a network of neighboring shapes
neighbors <- poly2nb(shapes) 

#create weights based off of a shapes neighbors
weights <-nb2listw(neighbors, style="B", zero.policy=TRUE)

#local moran's test for autocorrelation 
morans <- localmoran(shapes$counts, weights, zero.policy=TRUE)

#bind moran's pvalues back to the data
data$dbf$pval <- morans[,5]

#rewrite data
write.dbf(data, "counts.dbf")