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

Roger Bivand Roger.Bivand at nhh.no
Thu Jul 12 21:31:30 CEST 2012 

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

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