[R-sig-Geo] LocalG in spdep to calculate Gi*

[Babak Naimi]

Hi Juliann,

You can use lisa function in usdm package which I implemented for calculating some of local statistics including local Gi* for raster data.

Hope this would help.

Cheers,
Babak


----------------------------------------------------------------------

Message: 1
Date: Sun, 10 Mar 2013 11:48:00 -0700 (PDT)
From: Empty Empty <phytophthorasb at yahoo.com>
To: "r-sig-geo at r-project.org" <r-sig-geo at r-project.org>
Subject: [R-sig-Geo] LocalG in spdep to calculate Gi*
Message-ID:
        <1362941280.23258.YahooMailNeo at web142503.mail.bf1.yahoo.com>
Content-Type: text/plain

Hi,

I'm still a beginner at R and I'm trying also to learn some statistical techniques I haven't used before. Any suggestions will be appreciated.
???????? I have fairly large data sets (1 km population and climate data at a continent scale) that were originally raster. I couldn't figure out how to do local Gi* on a raster? -I'd love to know if there is a way to do it- so we've converted to points (and subset one country for learning how to do it). So my points are on a regular grid and I'm trying to use localG in spdep.
???????? What I have figured out is that I need to calculate neighbors (various options, e.g. knearneigh, but which one is appropriate depends on data - any guidance on determining which is appropriate?), then make that into an nb object (e.g. knn2nb), then make that into a spatial weights matrix or listw object (e.g.nb2listw, again pointers toward figuring out style would be welcome), then I can run localG. Am I on track so far?
??????? Then I want to somehow calculate statistical significance and plot the Gi* as a map. With so many points and comparisons, how do you figure out a critical value? Is it even possible/meaningful to calculate significance? The output R is the raw Gi* value (not a z-score), right?
??????? What I have below (gleaned from R help and some course notes I found online) worked on the subset and on the full continent data set and I got results that look reasonable and match up well with what my colleague did with another method. I guess that means my questions are more statistical than R, I apologize if this isn't the right place to ask those questions.
??????? I think that I could make the script more concise by not making each step into its own object, but I wasn't sure if I would need them again and each step takes a long time on the full data set. What is considered best practice in this regard?

library(spdep)
library(RANN)
library(rgdal)
library(classInt)

ethipt<-readOGR(".","ethipt")

coords<-coordinates(ethipt)
eth.knn <- knearneigh(coords, k=4,longlat = F, RANN= T) #calculate neighbors (still need to determine appropriate method & K)

eth.gal.nb<-knn2nb(eth.knn)? #make nb object

eth.listw <- nb2listw(include.self(eth.gal.nb), style="W")??? #give spatial weights to the neighbors list as listw object (still need to check appropriate "style")

ethG<-localG(ethipt$GRID_CODE,eth.listw)???? #calculate Gi*

#plot critical values of Gi*

sig3 <- classIntervals(ethG, n=3, style="fixed",
?????????????????????? fixedBreaks=c(min(ethG),-3.886, 3.886, max(ethG)))???? #3.886 isn't right CV -too many points? ?
cs<-c('blue','white','red') # blue = below -3, white = intermediate, red = above 3

sigcol <- findColours(sig3, cs)? # Map classes to colors

plot(ethipt, col=sigcol, axes=TRUE)
title("Local Gi* with 3.886 as critical value")
legend("left", fill=attr(sigcol, 'palette'), legend = names(attr(sigcol, 'table')))

Thanks a lot for any suggestions,

Juliann
        [[alternative HTML version deleted]]