[R-sig-Geo] question on raster Moran's I statistical significance

[Roger Bivand]

On Sat, 23 Oct 2021, Gabriel Cotlier wrote:

> Hello,
>
> I would like to estimate the pseudo p-value of Moran's I from raster data
> layer, by converting it to polygon data in order to be able to directly use
> the function  spdep::mc.moran() and getting in one step not only value of
> Moran'sI but its statistical significance without using other simulations
> but only that in the funcion spdep::mc.moran() as an alternative I have
> used the code below. However, since with the initial version of the code I
> got  several errors such as :
>
> Error in moran. mc(My_raster, lw.l, 999) : r is not a numeric vector
>
> So I converted "r" to a vector using as.vector() in :
>
> M <- moran.mc(as.vector(r), lw.l, 999)

First, if you have the spatial weights object lw.l, spdep::moran.test(), 
using analytical rather than bootstrap inference, is almost certainly more 
efficient, and yields almost always the same outcome with 
randomisation=TRUE.

>
> However then appeared another the error:
>
> Error in na.fail.default(x) : missing values in object
>

Please see the na.action= argument to moran.mc() and moran.test(); the 
default is na.fail, but can be na.exclude, which then also subsets the 
listw spatial weights object. You can also subset l before creating the 
neighbour object:

library(raster)
library(spdep)
r <- raster(nrows=10, ncols=10)
values(r) <- c(1:(ncell(r)-1), NA)
l <- rasterToPolygons(r, na.rm=TRUE) # note that this subsets by default
dim(l) # [1] 99  1
nb.l <- poly2nb(l, queen=FALSE)
nb.l # here no no-neighbour observations, no need for zero.policy=
lw.l <- nb2listw(nb.l, style="W")
moran.test(l$layer, lw.l) # default randomisation=TRUE



> Finally these errors disappeared with the following version of the code--
> see below -- the version of the code below runs OK--unless for the input
> raster I have used-- with no more errors or problems.
>

Because rasterToPolygons() removes NA observations by default.

Hope this clarifies,

Roger

> Therefore my question is : if instead of using arguments of the function
> spdep::mc.moran() such as na.action = na.omit or other --which I have tried
> and could not make work-- is it a valid way to solve the above errors  in
> the way as presented in the code below?
>
> I would appreciate it a lot if there is a possibility of knowing what could
> be wrong in the code below?
>
> I suspect it could be possible the problem is related to the way in which
> the raster  data layer is converted to polygon data  im the code to be
> input to the function spdep::mc.moran() , but will appreciate any guidance
> and assistance in the issue to get working in the correct manner.
>
> Thanks a lot.
> Kind regards,
>
> ################### Estimate of Moran's I and its p-value from a raster
> layer ##################################
> library(raster)
> library(spdep)
> library(easycsv)
>
> rm(list = ls())
>
> 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)
> dev.off()
> ################################################################################################
>
>

-- 
Roger Bivand
Emeritus Professor
Department of Economics, Norwegian School of Economics,
Postboks 3490 Ytre Sandviken, 5045 Bergen, Norway.
e-mail: Roger.Bivand using nhh.no
https://orcid.org/0000-0003-2392-6140
https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en