Is Moran’s I robust?

4 May, 2012 (12:15) | R, Spatial analysis | By: Manuel Gimond

Moran’s I (both global and local) is a measure of spatial autocorrelation. ¬†One form of its equation looks like this:

  I = \frac{N}{\sum_{i} \sum_{j} w_{ij}} \frac{\sum_{i} \sum_{j} w_{ij} (X_i - \bar X)(X_j - \bar X)}{\sum_{i} (X_i - \bar X )^2}

Note the residual and spread terms. As such, Moran’s I can be sensitive to outliers. The following example demonstrates this point:

# Load NC data
NC <- readShapePoly(system.file("etc/shapes/sids.shp", package="spdep")[1],
      ID="FIPSNO", proj4string=CRS("+proj=longlat +ellps=clrk66"))
rn <- sapply(slot(NC, "polygons"), function(x) slot(x, "ID"))
NB <-"etc/weights/", package="spdep")[1],
n  <- length(NB)
# Define some variables
x.norm <- rnorm(n)     # Normal distribution
rho    <- 0.3          # autoregressive parameter
W      <- nb2listw(NB) # Generate spatial weights
# Generate spatially autocorrelated datasets 
# (one normally distributed the other skewed) <- invIrW(W, rho) %*% x.norm # Generate autocorrelated values <- exp( # Transform orginal data to create a 'skewed' version

You can map the simulated data as follows:

# Append data to NC polygon object
auto_sk <- as.vector(
auto_no <- as.vector(
rn      <- as.integer(rn)
NCa     <- spCbind(NC,auto_sk)
NCa     <- spCbind(NCa,auto_no)
NCa     <- spCbind(NCa,rn)   # Needed if the GAL file is to be used in GeoDa
# Plot map (skewed data in this example)
brks <- quantile(auto_sk, seq(0,1,1/10))
cols <- grey.colors(length(brks)-1, 0.95, 0.55, 2.2)
plot(NCa, col=cols[findInterval(auto_sk, brks, all.inside=TRUE)])

# Compute Moran's I for both datasets
I.norm = moran.test(, listw=W)
I.skew = moran.test(, listw=W)
# The above approach makes assumptions about the Moran's I when in fact, it may not be.
# To correct for this, we can generate an empirical sampling distribution function
# using a Monte Carlo procedure.
MCI.norm <-, listw=W, nsim=9999)
MCI.skew <-, listw=W, nsim=9999)
# Compare the P-values

The p-values from the Monte Carlo simulations are 0.0133 and 0.1672 for the normal and skewed data respectively. The test results under different data distributions provide us with two distinct conclusions at a 5% significance level!

oopar <- par(mfrow=c(1,2), mar=c(4,4,3,2)+0.1)

NOTE: To export the spatial data to a shapefile for further analysis in GeoDa:

# Write map object to shapefile
writePolyShape(NCa,"NCa",factor2char = F)

–> If using Geoda, the weight file can be found under C:\Program Files\R\R-2.12.2\library\spdep\etc\weights