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

[Roger Bivand]

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

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