[R-sig-Geo] How to handle NA-values in raster-based Geary´s C test?

Kerstin Traut kerstin.traut at uni-jena.de
Thu Jul 12 22:16:55 CEST 2012 


-----Ursprüngliche Nachricht-----
Von: Roger Bivand [mailto:Roger.Bivand at nhh.no] 
Gesendet: Donnerstag, 12. Juli 2012 21:32
An: Kerstin Traut
Cc: r-sig-geo at r-project.org
Betreff: Re: [R-sig-Geo] How to handle NA-values in raster-based Geary´s C
test?

On Thu, 12 Jul 2012, Kerstin Traut wrote:

> Zitat von Roger Bivand <Roger.Bivand at nhh.no>:
>
>> On Thu, 12 Jul 2012, Kerstin Traut wrote:
>> 
>>> Zitat von Roger Bivand <Roger.Bivand at nhh.no>:
>>> 
>>>> On Thu, 12 Jul 2012, Kerstin Traut wrote:
>>>> 
>>>>> Zitat von Roger Bivand <Roger.Bivand at nhh.no>:
>>>>> 
>>>>>> On Thu, 12 Jul 2012, Kerstin Traut wrote:
>>>>>> 
>>>>>>> Zitat von Roger Bivand <Roger.Bivand at nhh.no>:
>>>>>>> 
>>>>>>>> On Thu, 12 Jul 2012, Roger Bivand wrote:
>>>>>>>> 
>>>>>>>>> On Thu, 12 Jul 2012, Kerstin Traut wrote:
>>>>>>>>> 
>>>>>>>>>> Zitat von Roger Bivand <Roger.Bivand at nhh.no>:
>>>>>>>>>> 
>>>>>>>>>>> On Thu, 12 Jul 2012, Kerstin Traut wrote:
>>>>>>>>>>> 
>>>>>>>>>>>> Hi,
>>>>>>>>>>>> 
>>>>>>>>>>>> I have a question on testing spatial autocorrelation on 
>>>>>>>>>>>> raster data including NA-values. In particular, I like to 
>>>>>>>>>>>> calculate Moran´s I and Geary´s C indices by using inverse 
>>>>>>>>>>>> distance weighting matrices.
>>>>>>>>>>>> Calculating Moran´s I with moran.test works fine, because 
>>>>>>>>>>>> the function contains the option "na.action=na.pass". 
>>>>>>>>>>>> Unfortunately, the function geary.test does not contain 
>>>>>>>>>>>> this option. The problem is that geary.test needs an nb 
>>>>>>>>>>>> argument for which I cannot perform for example 
>>>>>>>>>>>> nb<-ifelse(value, nb, NA), because it would change the
nb-format to a list-format.
>>>>>>>>>>>> Is there a way to ignore NA-values in geary.test? Or how 
>>>>>>>>>>>> can I set NA values in a neigbors (nb) list?
>>>>>>>>>>> 
>>>>>>>>>>> The simple answer is to subset the variable of interest and 
>>>>>>>>>>> the nb object to omit NA observations (which is what 
>>>>>>>>>>> moran.test does internally). Alternatively, you could 
>>>>>>>>>>> contribute a patch to geary.test(), if you like.
>>>>>>>>>>> 
>>>>>>>>>>> Hope this helps,
>>>>>>>>>>> 
>>>>>>>>>>> Roger
>>>>>>>>>> 
>>>>>>>>>> Thanks for your answer. Unfortunately, I think that 
>>>>>>>>>> subsetting is in my case not the right solution, because I 
>>>>>>>>>> need to define the amount of neighbors in cell2nb on the 
>>>>>>>>>> basis of a rectangular extent (e.g. 50 x 50 neighbors). 
>>>>>>>>>> Subsetting this extent by deleting NA-values would destroy the nb
list format required by geary.test.
>>>>>>>>>> Is there any solution for subsetting an nb list without 
>>>>>>>>>> destroying the nb format?
>>>>>>>>> 
>>>>>>>>> Yes, there are subset methods for nb and listw objects. To 
>>>>>>>>> ensure data integrity, use dnearneigh() instead of cell2nb(), 
>>>>>>>>> and use a SpatialPointsDataFrame to hold your data. In that 
>>>>>>>>> way, you make sure that the subset condition is applied to the
same observations.
>>>>>>>> 
>>>>>>>> The latter advice may not work if your coordinates are planar 
>>>>>>>> and you really mean to set torus=TRUE in cell2nb(). You'll have 
>>>>>>>> to get the data order right by deconstructing the region.id 
>>>>>>>> attribute of the nb object created by cell2nb() to ensure data
integrity.
>>>>>>>> 
>>>>>>>> Roger
>>>>>>> 
>>>>>>> I already tested dnearneigh(), but I like to avoid this function 
>>>>>>> because for large datasets cell2nb() calculates neighbors more 
>>>>>>> quickly.
>>>>>>> Thanks to your suggestion I was able to subset the nb list by 
>>>>>>> using nb.subset(). But I could not find a suitable function for 
>>>>>>> subsetting a listw. Using the subset() function results in the 
>>>>>>> following error
>>>>>>> message: "Not yet able to subset general weights lists". Do you 
>>>>>>> have an explaination?
>>>>>> 
>>>>>> Yes, because the changes may not be appropriate, depending on how 
>>>>>> the general weights were defined. In your case, subset the nb 
>>>>>> object and the list of idw weights, and rebuild the listw object 
>>>>>> with the subsetted objects.
>>>>>> 
>>>>>> Roger
>>>>> 
>>>>> That´s what I already tried, but the error message says "neighbors 
>>>>> and glist do not conform"??
>>>> 
>>>> Right, you'll need to rebuild from the subsetted nb object from the
>>>> nbdists() step (and subsetted coordinates), because the subsetted 
>>>> idw list will have weights for excluded observations.
>>>> 
>>>> Roger
>>> 
>>> Ok, I hope I understood everything right. You mean, firstly, I need 
>>> to subset my nb list and afterwards, I need to calculate the 
>>> nb-distances using nbdists()?
>> 
>> Yes, that will give consistent results when the listw object contains 
>> general weights. The same would apply to moran.test() with a listw 
>> object with general weights, I would have thought based on the code 
>> of subset.listw(), so I suggest not using na.pass in the moran.test() 
>> call, and doing the same reconstruction before calling the tests. I'd 
>> be interested in seeing an example of data (value, nb.idw, 
>> coordinates(subset)),  that run successfully through moran.test() 
>> because it should also fail on general weights.
>> 
>> Roger
>
> Finally, after some more failed trials the geary.test() is working for 
> my dataset. Thank you so much for your suggestions!
> I tested moran.test() with the consistent method against the na.pass
method. 
> For the first method I got an I-value of 0.7952, for the second a 
> Moran´s I value of 0.79477. Do you know, how to explain this small 
> (but existing) difference?

Using na.action=na.pass allows NA values, and will give a different result
at least, I think, because na.pass does not adjust n (subtracting NA
observations), nor does it subset. If you use na.action=na.omit, n will be
reduced as in the consistent method, but subset.listw() fails because the
"mode" attribute of the weights component of the listw object is not
"binary". I think here that you should rely on subsetting manually.

Hope this clarifies,

Roger

Thank you very very much! Now everything is clear!

Kind regards,
Kerstin

>
> Kerstin
>
>>> 
>>> Kerstin
>>> 
>>> 
>>>>> 
>>>>> Kerstin
>>>>> 
>>>>>>> 
>>>>>>> Kerstin
>>>>>>> 
>>>>>>> 
>>>>>>>>> 
>>>>>>>>> Roger
>>>>>>>>> 
>>>>>>>>>> 
>>>>>>>>>> Kind regards,
>>>>>>>>>> Kerstin
>>>>>>>>>> 
>>>>>>>>>> 
>>>>>>>>>>>> Thank you very much,
>>>>>>>>>>>> Kerstin
>>>>>>>>>>>> 
>>>>>>>>>>>> ---
>>>>>>>>>>>> 
>>>>>>>>>>>> file<-raster(test.tif")
>>>>>>>>>>>> mask<-extent(file, 1, 50, 1, 50) 
>>>>>>>>>>>> subset<-crop(file,mask,test.tif", overwrite=TRUE) value <- 
>>>>>>>>>>>> as.matrix(subset)
>>>>>>>>>>>> 
>>>>>>>>>>>> # Creation of a list of integer vectors giving the region 
>>>>>>>>>>>> id numbers for the neighbors within the grid extent nb <- 
>>>>>>>>>>>> cell2nb(50, 50, type="queen", torus=TRUE)
>>>>>>>>>>>> 
>>>>>>>>>>>> # Calculation of the distances along the links in the 
>>>>>>>>>>>> neighbous list dist <- nbdists(nb, coordinates(subset), 
>>>>>>>>>>>> longlat=TRUE)
>>>>>>>>>>>> 
>>>>>>>>>>>> # Converting the distances to inverse distance values 
>>>>>>>>>>>> inv.dist <- lapply(dist, function (x) 1/(x*100))
>>>>>>>>>>>> 
>>>>>>>>>>>> # Supplementing the neighbors list with spatial weights for 
>>>>>>>>>>>> a coding scheme nb.idw <- nb2listw(nb, glist=inv.dist, 
>>>>>>>>>>>> style="W")
>>>>>>>>>>>> summary(unlist(nb.idw$weights))
>>>>>>>>>>>> 
>>>>>>>>>>>> # Moran's test (two-sided)
>>>>>>>>>>>> moran.global <- moran.test(value, listw=nb.idw, zero.policy 
>>>>>>>>>>>> = TRUE, na.action=na.pass, alternative="two.sided") 
>>>>>>>>>>>> moran.global
>>>>>>>>>>>> 
>>>>>>>>>>>> # Geary's test (two-sided)
>>>>>>>>>>>> geary.global <- geary.test(value, listw=nb.idw, zero.policy 
>>>>>>>>>>>> = TRUE , alternative="two.sided") geary.global
>>>>>>>>>>>> 
>>>>>>>>>>>> -----------------------------------------------------------
>>>>>>>>>>>> ----- This message was sent through 
>>>>>>>>>>>> https://webmail.uni-jena.de
>>>>>>>>>>>> 
>>>>>>>>>>>> _______________________________________________
>>>>>>>>>>>> R-sig-Geo mailing list
>>>>>>>>>>>> R-sig-Geo at r-project.org
>>>>>>>>>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>>>>>>>>> 
>>>>>>>>>>> --
>>>>>>>>>>> Roger Bivand
>>>>>>>>>>> Department of Economics, NHH Norwegian School of Economics, 
>>>>>>>>>>> Helleveien 30, N-5045 Bergen, Norway.
>>>>>>>>>>> voice: +47 55 95 93 55; fax +47 55 95 95 43
>>>>>>>>>>> e-mail: Roger.Bivand at nhh.no
>>>>>>>>>>> 
>>>>>>>>>> 
>>>>>>>>>> 
>>>>>>>>>> 
>>>>>>>>>> -------------------------------------------------------------
>>>>>>>>>> --- This message was sent through https://webmail.uni-jena.de
>>>>>>>>> 
>>>>>>>>> 
>>>>>>>> 
>>>>>>>> --
>>>>>>>> Roger Bivand
>>>>>>>> Department of Economics, NHH Norwegian School of Economics, 
>>>>>>>> Helleveien 30, N-5045 Bergen, Norway.
>>>>>>>> voice: +47 55 95 93 55; fax +47 55 95 95 43
>>>>>>>> e-mail: Roger.Bivand at nhh.no
>>>>>>>> 
>>>>>>> 
>>>>>>> 
>>>>>>> 
>>>>>>> ----------------------------------------------------------------
>>>>>>> This message was sent through https://webmail.uni-jena.de
>>>>>> 
>>>>>> --
>>>>>> Roger Bivand
>>>>>> Department of Economics, NHH Norwegian School of Economics, 
>>>>>> Helleveien 30, N-5045 Bergen, Norway.
>>>>>> voice: +47 55 95 93 55; fax +47 55 95 95 43
>>>>>> e-mail: Roger.Bivand at nhh.no
>>>>>> 
>>>>> 
>>>>> 
>>>>> 
>>>>> ----------------------------------------------------------------
>>>>> This message was sent through https://webmail.uni-jena.de
>>>> 
>>>> --
>>>> Roger Bivand
>>>> Department of Economics, NHH Norwegian School of Economics, 
>>>> Helleveien 30, N-5045 Bergen, Norway.
>>>> voice: +47 55 95 93 55; fax +47 55 95 95 43
>>>> e-mail: Roger.Bivand at nhh.no
>>>> 
>>> 
>>> 
>>> 
>>> ----------------------------------------------------------------
>>> This message was sent through https://webmail.uni-jena.de
>> 
>> --
>> Roger Bivand
>> Department of Economics, NHH Norwegian School of Economics, 
>> Helleveien 30, N-5045 Bergen, Norway.
>> voice: +47 55 95 93 55; fax +47 55 95 95 43
>> e-mail: Roger.Bivand at nhh.no
>> 
>
>
>
> ----------------------------------------------------------------
> This message was sent through https://webmail.uni-jena.de

--
Roger Bivand
Department of Economics, NHH Norwegian School of Economics, Helleveien 30,
N-5045 Bergen, Norway.
voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand at nhh.no

More information about the R-sig-Geo mailing list