[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 18:04:11 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, 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



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

More information about the R-sig-Geo mailing list