[R-sig-Geo] question on raster Moran's I statistical significance
Roger Bivand
Roger@B|v@nd @end|ng |rom nhh@no
Mon Oct 25 13:43:45 CEST 2021
On Mon, 25 Oct 2021, Roger Bivand 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?
An afterthought - are there any diagonal linear patterns? You have chosen
rook not queen, so a diagonal line of cells remaining after dropping
NAs will generate singleton observations along that line,even though they
would have been queen neighbours as they touch at one point. What is the
coordinnate reference system of the input raster (crs(r))?
Roger
>
> 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
More information about the R-sig-Geo
mailing list