[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 17:42:13 CEST 2012


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?

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



More information about the R-sig-Geo mailing list