Raster – Moran’s I Statistical Significance with spdep mc.moran Function

autocorrelationrasterspdep

I would like to obtain Moran's I estimation together with its p-value from a raster data layer –which has na value– including a figure showing the simulation distribution by mean of spdep::mc.moran() function.

I would like to know whether the code bellow implies the right manner to convert the rater to polygon as required to be input to the function spdep::mc.moran(), as well as whether the code bellow is an appropriate manner to alternatively solve the problem ?

library(raster)
library(spdep)

r<- raster(file.choose())
l <- rasterToPolygons(r)

nb.l <- poly2nb(l, queen=FALSE)
lw.l <- nb2listw(nb.l, style="W",zero.policy=TRUE)

v<-r[!is.na(r)]
M <- moran.mc((v), lw.l, 999, zero.policy=TRUE)
M
  
plot(M, main=sprintf("Original, Moran's I (p = %0.3f)", M$p.value), xlab="Monte-Carlo Simulation of Moran I", sub =NA, cex.lab=1.5, cex.axis=1.5, cex.main=1.5, cex.sub=1.5)
abline(v=M$statistic, col="red", lw=2)

Best Answer

This looks perfectly correct to me - you can reproduce the value of Moran(r) if you use queen=TRUE in the neighbour calculation and use style="B" in the weighted neighbour calculation, so that all the weights are 1. With style="W" the sum of weights is 1 for each feature.

By default raster::Moran uses a weight matrix:

 w = matrix(c(1, 1, 1, 1, 0, 1, 1, 1, 1), 3, 3)

which corresponds to queen's adjacency and a binary adjacency. Using my standard r = raster(matrix(1:12, 4,3)) quick test raster, and following your code:

> nb.l <- poly2nb(l, queen=TRUE)
> lw.l <- nb2listw(nb.l, style="B",zero.policy=TRUE)
> v<-r[!is.na(r)]
> M <- moran.mc((v), lw.l, 999, zero.policy=TRUE)
> M

    Monte-Carlo simulation of Moran I

data:  (v) 
weights: lw.l  
number of simulations + 1: 1000 

statistic = 0.33205, observed rank = 994, p-value = 0.006
alternative hypothesis: greater

> Moran(r)
[1] 0.3320473
Related Question