The Moran statistic

An overall measure of spatial autocorrelation

The Moran statistic provides a one-number overall measure of spatial autocorrelation.  The statistic is broadly analogous to the ordinary correlation coefficient, with a numerator that measures the extent to which adjacent points (as described by the elements of W) have similar deviations about the mean of the data, and a denominator that standardizes that quantity to reflect the scale or variability of the variable being examined.

The Moran statistic

Start by mapping a couple of gridded climate variables to illustrate the overall trends:

library(maptools
library(RColorBrewer)
library(spdep)

# read a county outline shape file
orotl.shp <- readShapeLines(file.choose(),
     proj4string=CRS("+proj=longlat"))
# read and attach some gridded climate data to illustrate overall pattern
orgrid <- read.csv(file.choose())

attach(orgrid)

# set up for maps of individual variables
brks <- c(0,10,20,50,100,200,500,1000,9000)
plotvar <- pjan # pick a variable to plot
plottitle <- "January Precipitation"
nclr <- 8 # number of colors
plotclr <- brewer.pal(nclr,"BuPu") # get the colors
#plotclr <- plotclr[nclr:1] # reorder colors if appropriate
colornum <- as.integer(cut(plotvar, brks, include.lowest=T))
colcode <- plotclr[colornum] # assign color

# plot the county outline shape file and the variable
plot(orotl.shp, xlim=c(-125, -115), ylim=c(42,47), axes=T)
points(lon, lat, pch=15, col=colcode, cex=1.2)
text(-120, 46.5, plottitle)

# add legend
legend(-116.15, 44, legend=leglabs(brks), fill=plotclr, cex=0.6, bty="o")
plot(orotl.shp, add=T)

detach(orgrid)

Now attach the station dataset, and get the distance-based neighbors

attach(orstationc)

# station location matrix
orstationc.loc <- cbind(lon, lat)

# plot the county outline shape file and the variable
plot(orotl.shp, xlim=c(-125, -115), ylim=c(42,47), axes=T)

# add climate station points
points(orstationc.loc, pch=16, cex=1.2)

# distance neighbors
d <- 80
orstationc.neighbors.dist <- dnearneigh(orstationc.loc, 0, d, longlat=TRUE)

plot(orstationc.neighbors.dist, orstationc.loc, add=T, lwd=2, col="blue")
text(-120, 46.5, "Distance-based neighbors (d = 80 km)")

Plot one of the station dataset variables:

brks <- c(0,10,20,50,100,200,500,1000,9000)
plotvar <- pjan # pick a variable to plot
plottitle <- "January Precipitation"
nclr <- 8 # number of colors
plotclr <- brewer.pal(nclr,"BuPu") # get the colors
#plotclr <- plotclr[nclr:1] # reorder colors if appropriate
colornum <- as.integer(cut(plotvar, brks, include.lowest=T))
colcode <- plotclr[colornum] # assign color

# plot the county outline shape file and the variable
plot(orotl.shp, xlim=c(-125, -115), ylim=c(42,47), axes=T)
points(lon, lat, pch=16, col=colcode, cex=2)
points(orstationc$lon, orstationc$lat, cex=2)
text(-120, 46.5, plottitle)

# add legend
legend(-116.15, 44, legend=leglabs(brks), fill=plotclr, cex=0.6, bty="o")
plot(orotl.shp, add=T)

Finally, get the Moran statistic:

# Moran statistic
moran.test (log10(pjan),
    nb2listw(orstationc.neighbors.dist, style="W"))

The two arguments of the moran.test() function are the name of the variable being analyzed, and the results of applying the nb2listw() function to a spatial-neighbors object, which formats it in an appropriate fashion.  The style argument of nb2listw() function indicates how to treat the individual spatial weights (e.g. either as basic binary coding (B) or as row-standardized values (W)).

An interesting plot can be constructed that shows the relationship between  the values of the variable being analyzed plotted relative to the spatially lagged values of that variable, as represented by the spatial neighbors.  This plot can examined in much the same way as a typical scatter diagram, including adding a straight-line summary, and analyzing the influence of individual observations on that relationship.

moran.plot(pjul, nb2listw(orstationc.neighbors.dist),
    labels=as.character(station), pch=19)