I would like to estimate the Moran's I coefficient for raster data together with the statistical significance of the spatial autocorrelation obtained.

I found that the raster package function Moran() although calculates the spatial autocorrelation index it apparently does not give directly the statistical significance of the results obtained: https://search.r-project.org/CRAN/refmans/raster/html/autocor.html

Could it be possible to obtain the statistical significance of the results with either a raster package or similar one?

1 Answer 1


You could do a monte-carlo test of I>0:

First lets create a very correlated raster:

> r = raster(matrix(1:(50*50),50,50))
> Moran(r)
[1] 0.9694908

And now do 99 Moran's I of rasters that are random samplings of those values:

> M99 = sapply(1:99, function(i){v = r; v[]=sample(r[]);Moran(v)})

And let's see the distribution:

> hist(M99)
> range(M99)
[1] -0.02713834  0.02061854

And its clear that the Moran(r) is way outside the range of the simulations, so we can reject the null that the data are uncorrelated with respect to random sampling/rearrangement of the data.

To get the approximate pseudo p-value, see where the Moran stat for the data ranks amongst the simulations. Suppose the Moran stat for the data was 0.018 (using the 0.969 example from my code is a bit extreme), then compute the rank:

> rank(c(0.018, M99))[1]
[1] 98

Which shoes that 0.018 ranks 98 out of the 100 values (99 sims + itself). Hence reject H0 (no spatial autocorrelation under random arrangement hypothesis) with approximate p = 0.98.

If you do more simulations, then the general case is:

> R/(length(M99)+1)
[1] 0.98

for any number of simulations in M99.

This is what moran.mc from the spdep package does with data for polygons or other general neighbourhood structures.

Alternatively convert your raster to a grid of polygons and use the spdep functions with either 4-way or 8-way neighbours (or beyond...)

  • Thanks a lot for the answer is indeed very helpful ! Just a question the value of 99 is valid for any raster layer size or should be related to the size of the raster layer?
    – Gab
    Commented Oct 20, 2021 at 17:20
  • That's the number of monte-carlo simulations to do. The more you do, the more precise is your estimate of the (pseudo) p-value. 99 simulations plus your observational data mean your p-value is the rank of the observation out of 100 - in this case 100/100, and I can say p > 0.99. If I did 999 sims and it was still largest, that would be p > 0.999. If I did 999 sims and it was 994th, then p>0.994. See ?moran.mc in spdep and general MC test theory
    – Spacedman
    Commented Oct 21, 2021 at 7:47
  • Thanks a lot for the clear explanation.
    – Gab
    Commented Oct 21, 2021 at 21:37
  • 1
    is there any chance to get for the same example code another line for obtaining the correspondent p-value to that example?
    – Gab
    Commented Oct 22, 2021 at 13:42
  • Just a question to clarify myself with regards to pseudo p-value estimation: in my case Moran(r) = 0.6461546, so the pseudo p-value would be 1, right? since rank(c(0.6461546, M99))[1] is 100, so pseudo would p-value = 1 then, H0 is rejected, right ?
    – Gab
    Commented Oct 22, 2021 at 14:36

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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