R Moran Index – Statistical Significance of Moran I

Question

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?

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...)