[R-sig-Geo] Standardizing Moran's I

[Giancarlo Sadoti]

Greeting list members,

I'm analyzing the global spatial autocorrelation of a continuous vegetation attribute (MSH) among multiple (15-30) plots (secondary sampling units) across 76 sites (primary sampling units).  I've used the Moran's I permutation test (moran.mc) in the spdep package to do this.  I've included my code from one of my 76 sites:

----------
> library(spdep)
> head(data)
    EAST   NORTH MSH
1 274276 4818867  70
2 274300 4818942  70
3 274249 4818900  60
4 274258 4818900  60
5 274299 4818900  66
6 274304 4818900  66

> dist=as.matrix(dist(cbind(data$EAST,data$NORTH)))
> inv.dist=1/dist
> diag(inv.dist)=0
> lw=mat2listw(inv.dist)
> m=moran.mc(data$MSH,lw,nsim=1000)
> m

        Monte-Carlo simulation of Moran's I

data:  data$MSH 
weights: lw  
number of simulations + 1: 1001 
 
statistic = 0.1704, observed rank = 999, p-value = 0.001998
alternative hypothesis: greater 

> str(m)
List of 7
 $ statistic  : Named num 0.17
  ..- attr(*, "names")= chr "statistic"
 $ parameter  : Named num 999
  ..- attr(*, "names")= chr "observed rank"
 $ p.value    : num 0.002
 $ alternative: chr "greater"
 $ method     : chr "Monte-Carlo simulation of Moran's I"
 $ data.name  : chr "data$MSH \nweights: lw  \nnumber of simulations + 1: 1001 \n"
 $ res        : num [1:1001] 0.0193 -0.00302 -0.00983 -0.00605 -0.03424 ...
 - attr(*, "class")= chr [1:2] "htest" "mc.sim"

------------

I'm investigating whether the spatial autocorrelation of the attribute across the 76 sites is a strong predictor of an ecological response at each site.  Because sample sizes (and expected values of I) vary across sites (-1/[N-1]; N varies between 15-30), the observed Moran's I at each site appears to be an inappropriate statistic for comparison across sites.   I originally proposed the standardization of Moran's I to z-values using the observed Moran's I (I) and the mean (MEAN_I) and sd of I (SD_I) of the 1,000 simulations:

z = (I - MEAN_I) / SD_I

However, as suspected by a reviewer (and confirmed by looking at my data), most (if not all) of the 76 site distributions of the 1,000 expected values of I deviate from the assumption of a normal distribution (tested using Shapiro-Wilks and Kolmogorov-Smirnov tests); most are right-skewed.  Likewise, the p-value of the statistic is an inappropriate standardization due to both positive and negative autocorrelation across sites.  I've considered using the rank of the statistic within the distribution of expected values, but this seems problematic in a way I can't quite articulate.

Can anyone offer any suggestions for an alternative way to standardize Moran's I values in a situation such as this?  Your input will be greatly appreciated.  My apologies in advance if this has been discussed previously.

Giancarlo