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

Mon Oct 25 14:27:01 CEST 2021

Thanks a lot.
I will try spdep::n.comp.nb()  and read the link you suggested to me
on "Proximity and aerial data" which seems very interesting.
Maybe could it be due to the fact that the data layer of the city is
traversed by a river, which I set to na since it was not relevant for the
analysis, thus maybe from there the message "reginos with no links"?
However when I do :
> nb.l
Neighbour list object:
Number of regions: 636410
Number of nonzero links: 5060866
Percentage nonzero weights: 0.001249542
Average number of links: 7.95221
8 regions with no links:
103556 104030 104502 104970 170752 290616 291628 520687

Then, even with regions with no links, I have no error when setting
zero.policy = TRUE in nb2listw()

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

and then I can run moran.test() which gives some results :

> moran.test(l$dta, lw.l,zero.policy = TRUE) # default randomisation=TRUE

Moran I test under randomisation

data:  l$dta
weights: lw.l  n reduced by no-neighbour observations

Moran I statistic standard deviate = 1559.1, p-value <
2.2e-16
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance
     9.815120e-01     -1.571336e-06      3.963186e-07

Why could it be running apparently ok in this way ?

Thanks a lot again for your help and guidance.
Kind regards,
Gabriel Cotlier

On Mon, Oct 25, 2021 at 2:12 PM Roger Bivand <Roger.Bivand using nhh.no> wrote:

> On Mon, 25 Oct 2021, Gabriel Cotlier wrote:
>
> > Hello,
> >
> > I have tried the code for most of the dataset but only few make the
> > same error:
> >> nb.l
> > Neighbour list object:
> > Number of regions: 975434
> > Number of nonzero links: 3869290
> > Percentage nonzero weights: 0.0004066638
> > Average number of links: 3.966737
> > 79 regions with no links:
> > 3206 34632 47267 65155 81013 90934 91361 91427 98318
> > 98765 100054 100472 100884 100896 102508 102514 105207
> > 105937 106630 106965 114449 115479 122329 123614 126880
> > 128631 314459 350405 381283 397379 401275 412489 420698
> > 426384 428632 429752 432058 433230 440631 451195 451196
> > 452270 464327 467009 503242 504459 507625 510412 517865
> > 517916 524413 526766 535269 537320 537712 544576 549519
> > 551629 553413 565169 566421 567703 567791 568587 579209
> > 585464 625809 642179 648122 651096 652581 654065 708908
> > 777500 845322 846238 847153 848065 866391
> >> lw.l <- nb2listw(nb.l, style="W")
> > Error in nb2listw(nb.l, style = "W") : Empty neighbour sets found
> >
> > Could it be something related to the input data layer?
>
> Yes. Please also look at the output of spdep::n.comp.nb() if it works for
> a dataset this large, see ?n.comp.nb and examine the nc component (number
> of subgraphs). The singleton graphs have no neighbours, but there may be
> other subgraphs, you find their sizes by table(table(<comp.id
> component>)). This may well have come from the na.rm=TRUE in
> rasterToPolygons(), but is unavoidable. Maybe see also:
> https://keen-swartz-3146c4.netlify.app/area.html#
> and expand all code (button at top right) and search for n.comp.nc.
>
> Since you have singleton observations, zero.policy=TRUE is your only
> option unless you further subset by omitting all the observations for
> which card(nb.l) == 0L - since they have no neighbours, they do not inform
> the spatial process.
>
> Roger
>
> > How could it be possible to avoid this error for running the code?
> >
> > Thanks a lot for your guidance and help.
> >
> > Kind regards,
> > Gabriel Cotlier
> >
> >
> > On Mon, Oct 25, 2021 at 12:21 PM Gabriel Cotlier <gabiklm01 using gmail.com>
> > wrote:
> >
> >> Thank you very much for your response.
> >> Yes indeed it clarifies and helps a lot.
> >> Kind regards,
> >> Gabriel Cotlier
> >>
> >> On Mon, Oct 25, 2021 at 12:00 PM Roger Bivand <Roger.Bivand using nhh.no>
> wrote:
> >>
> >>> 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
> >>>
> >>
> >>
> >> --
> >> Gabriel Cotlier, PhD
> >> Haifa Research Center for Theoretical Physics and Astrophysics (HCTPA)
> >> University of Haifa
> >> Israel
> >>
> >>
> >>
> >
>
> --
> 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
>

	[[alternative HTML version deleted]]