2

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)

1 Answer 1

3

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
1
  • Thank you very much. Results now improved a lot.
    – Gab
    Oct 23, 2021 at 11:09

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.