[STATSGRASS] Re: [GRASS5] r.moran

[Roger Bivand]

On Fri, 24 Mar 2006, Wolf Bergenheim wrote:

> On Thu, 23 Mar 2006, Roger Bivand wrote:
> 
> >>
> >>>
> >>> The next bits are simple:
> >>>
> >>> my_nb <- poly2nb(my_shape)
> >>
> >> R says: Error: couldn't find function "poly2nb"
> >>
> >> perhaps something is missing?
> >>
> >> Ahh... I found it. help.search("poly2nb") talks about spdep.
> >> so I had to:
> >>
> >> library(spdep)
> >
> > No, I didn't do that on purpose just to see how you did! Very typically
> > users and developers work in an edit/run/edit/run/... mode, I tend to use
> > Vim to edit and source() to run (or copy&paste), very many power users use
> > Emacs Speaks Statistics which has a run buffer, and the same facilities
> > are available in RWinEdt and other editors. So a lot of things get
> > forgotten first time through ...
> >
> >>
> >> This is simply just great. I really wish I'd have had more time to study
> >> R. I guess now is as good a time as any.
> >
> > Being driven by pressing need is always good motivation, IMHO!
> >
> 
> (: Hmm is there any way to zoom in on the plot?

Use plot() again with the same arguments but add xlim= and ylim=, there is 
no interactive zoom on the standard output devices. You can grab the 
coordinates with locator().

> 
> >>
> >>>
> >>> maybe setting queen=FALSE to only choose neighbours sharing more than
> >>> just
> >>> a corner point.
> >>
> >> Hmm it complains that my_shape isn't a polygon list... ):
> >> What is the correct spell?
> >
> > Another slip on my part - I should have cast to a SpatialPolygons object:
> >
> > my_nb <- poly2nb(as(my_shape, "SpatialPolygons"))
> >
> > which will also run for a fair while.
> 
> Ah. ok, now it does indeed run!
> 
> >
> >>
> >>>
> >>> print(my_nb)
> >>> plot(my_nb, coordinates(my_shape))
> >>>
> >>> display the neighbours.
> >>>
> >>> From there it is just
> >>>
> >>> moran.test(my_shape$variable_of_interest, listw=nb2listw(my_nb, style="W",
> >>> zero.policy=TRUE), zero.policy=TRUE)
> >>>
> >>> roughly, depending on whether you have cells that have zero neighbours.
> >>> But we can get back to this once the shapefile is on board.
> 
> Yes indeed.  now I do get Moran's I and it looks just about right. I did 
> the same analysis in ArcMap, and the results look similar. The I value 
> differs a bit though... I wonder how critical this is.
> 
> In ArcMap I got Moran's I as 0.35 and with R I get 0.54. ArcMap gave a Z 
> value of 217.3 and R gave 52.2.

That is useful. I've got another conversation going off-list on the R side 
about the Arc Moran's I, there are issues with how the spatial weights are 
specified. We've managed to reproduce some values, but not all. It will be 
documented on the R-sig-geo archive (which is searchable through 
RSiteSearch() from the R prompt) before long, but if you could add your 
case, that would help a lot. 

What we need is the R commands used to create the weights plus the results
of running moran.test() on the data, and the equivalent (copy&paste from
the ArcGIS Python results popup window (the one with the progress bar at
the top with a good deal of meaningless text and number that are not
properly rounded)).

So far we know that for a points shapefile, the default in ArcGIS is to
join all points to all and use inverse distance, not standardising the
weights. Also the flat distance threshold can be reproduced. We haven't
tried a polygon layer yet.


If you could contribute matched argument sets with the same output, that 
would be great, but any documentation of what is going on will be 
valuable, either to show that spdep/GeoDa and ArcGIS give the same results 
or otherwise.

Thanks for your patience!

Roger

> 
> >>
> >> OK. Now I have the shapefile locked and loaded, now to just pull the
> >> moran numbers out of it (: and I will be solid.
> >>
> >
> > moran.test() returns an object of class "htest", which has a print method.
> > So what you see is the output of that. How do you need the output? The
> >
> > res <- moran.test(my_shape$variable_of_interest, listw=nb2listw(my_nb,
> >  style="W", zero.policy=TRUE), zero.policy=TRUE)
> > res # just to have a look
> > sink("tmp/output_file.txt")
> > res
> > sink()
> 
> Nice!
> 
> >
> > makes a text file copy. You can get it as HTML formatted too:
> >
> > library(R2HTML)
> > dir.create("tmp/morantest")
> > HTMLStart("tmp/morantest")
> > HTML(res)
> > HTMLStop()
> >
> 
> cool!
> 
> > but you'll need to edit it a bit afterwards (or copy & paste from a
> > browser to say oowriter). There is pretty good LaTeX support too, as
> > you'd expect with statisticians.
> >
> 
> That can come in handy as well.
> 
> > By the way, if names(my_shape) says for example:
> >
> > "children", "totalPop"
> >
> > and you want the percentage,
> >
> > my_shape$percent_children <- (100*my_shape$children)/my_shape$totalPop
> >
> > does the stuff.
> 
> Good to know.
> 
> >
> > Or if you want the local Moran's I in the same object, do str() on the
> > output object of localmoran() (a matrix, nb. different access syntax than
> > data frames) to work our which column you want - probably Z.Ii, put it in
> > my_shape, and do:
> >
> > writePolyShape(my_shape, "new_my_shape")
> >
> > to write out to a new shapefile. At this point simple working out the
> > appropriate row ID and writing out using functions in the RSQLite package
> > might be the most natural for you?
> >
> 
> I think I might go via the shape file. it seems more straight forward. (; 
> but I'll definitely look into the RSQLite package.
> 
> > One of the ideas Markus and Helena mention in their book is writing small
> > shell scripts from GRASS passing some values to R through environmental
> > variables to be used in an R script with the R interface hidden from the
> > user, because R can be run with a script from the command line.
> >
> 
> Heh, I expected as much (;
> 
> > So "get the local moran values for my polygons" could be run from the
> > GRASS prompt if enough of the options and other user decisions can be
> > handled sensibly. R can write to DBF too, and most external databases, but
> > the very range of choices makes writing a robust script hard. In the
> > polygon contiguity neighbour definition, we can actually pass out a
> > neighbour structure directly from GRASS, and R can read the GRASS
> > attribute databases directly, so it could be made quite "tight".
> >
> > Just thinking aloud ...
> 
> Hmm good ideas. (:
> 
> --Wolf
> 
> 

-- 
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, 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