Skip to content

How to handle NA-values in raster-based Geary´s C test?

15 messages · Roger Bivand, Kerstin Traut, djrwicks

Original
Kerstin Traut
# 
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?

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
Roger Bivand
#
On Thu, 12 Jul 2012, Kerstin Traut wrote:

            

      
      
      
    

        
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

            

      
      
      
    

        

      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Kerstin Traut
    

    

      
      #
    

  


  

      
Zitat von Roger Bivand <Roger.Bivand at nhh.no>:

            

      
      
      
    

        
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?

Kind regards,
Kerstin

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        
----------------------------------------------------------------
This message was sent through https://webmail.uni-jena.de

  
  


  


      

    

    

      
      

        


  

        

      Roger Bivand
    

    

      
      #
    

  


  

      
On Thu, 12 Jul 2012, Kerstin Traut wrote:

        

            

      
      
      
    

        
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.

Roger

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Roger Bivand
    

    

      
      #
    

  


  

      
On Thu, 12 Jul 2012, Roger Bivand wrote:

        

            

      
      
      
    

        
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

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Kerstin Traut
    

    

      
      #
    

  


  

      
Zitat von Roger Bivand <Roger.Bivand at nhh.no>:

            

      
      
      
    

        
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

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        
----------------------------------------------------------------
This message was sent through https://webmail.uni-jena.de

  
  


  


      

    

    

      
      

        


  

        

      Roger Bivand
    

    

      
      #
    

  


  

      
On Thu, 12 Jul 2012, Kerstin Traut wrote:

        

            

      
      
      
    

        
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 Traut
    

    

      
      #
    

  


  

      
Zitat von Roger Bivand <Roger.Bivand at nhh.no>:

            

      
      
      
    

        
That?s what I already tried, but the error message says "neighbors and  
glist do not conform"??

Kerstin

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        
----------------------------------------------------------------
This message was sent through https://webmail.uni-jena.de

  
  


  


      

    

    

      
      

        


  

        

      Roger Bivand
    

    

      
      #
    

  


  

      
On Thu, 12 Jul 2012, Kerstin Traut wrote:

        

            

      
      
      
    

        
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

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Kerstin Traut
    

    

      
      #
    

  


  

      
Zitat von Roger Bivand <Roger.Bivand at nhh.no>:

            

      
      
      
    

        
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()?

Kerstin

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        
----------------------------------------------------------------
This message was sent through https://webmail.uni-jena.de

  
  


  


      

    

    

      
      

        


  

        

      Roger Bivand
    

    

      
      #
    

  


  

      
On Thu, 12 Jul 2012, Kerstin Traut wrote:

        

            

      
      
      
    

        
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 Traut
    

    

      
      #
    

  


  

      
Zitat von Roger Bivand <Roger.Bivand at nhh.no>:

            

      
      
      
    

        
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?

Kerstin

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        
----------------------------------------------------------------
This message was sent through https://webmail.uni-jena.de

  
  


  


      

    

    

      
      

        


  

        

      djrwicks
    

    

      
      #
    

  


  

      
Hi Kerstin,

I have been suffering from a similar problem to you with NA values in my
raster dataset.  I was wondering if you might post or email me your revised
script so I can use it to troubleshooting along with the suggestions in this
thread?

Many thanks,

Dan.

--
View this message in context: http://r-sig-geo.2731867.n2.nabble.com/How-to-handle-NA-values-in-raster-based-Geary-s-C-test-tp7580439p7580454.html
Sent from the R-sig-geo mailing list archive at Nabble.com.

  
  


  


      

    

    

      
      

        


  

        

      Roger Bivand
    

    

      
      #
    

  


  

      
On Thu, 12 Jul 2012, Kerstin Traut wrote:

        

            

      
      
      
    

        
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 Traut
    

    

      
      #
    

  


  

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

        

            

      
      
      
    

        
nb-format to a list-format.

            

      
      
      
    

        
list format required by geary.test.

            

      
      
      
    

        
same observations.

            

      
      
      
    

        
integrity.

            

      
      
      
    

        
method.

            

      
      
      
    

        
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

            

      
      
      
    

        

      
      
    

            

      
      
      
    

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