[R-sig-Geo] Different Moran's I for same model residuals in different runs using moran.mc{spdep}

[Xingli Giam]

Thanks Roger for the reply and advice.

Rest assured I am taking my time with this. I have been spending the last week
looking at all the relevant examples, and running through the code given in the
helpfiles, as well as your spdep vignette. I am also using my own data as
examples to learn more about the different functions.

I realised my code for dnearneigh() is wrong! I used as.matrix instead of cbind
to create the matrix of long-lat coordinates!

That resulted in my weight lists object being inconsistent! I used cbind, and
different sessions returned the same Moran's I statistic (used lm.morantest on
your advice) on my same data and weights!

Thank goodness for the list! I am learning alot everyday!

Many thanks,
Xingli
(P.S. Thanks for the correlog{pgirmess} tip!)


Quoting Roger Bivand <Roger.Bivand at nhh.no>:

> On Sat, 31 Jan 2009, Xingli Giam wrote:
>
> > Hello List,
> >
> > I am using moran.mc{spdep} to extract values of Moran's I for a simple GLM
> > model residuals of binned distances.
> >
> > However, everytime I run moran.mc() in a new session, Moran's I values
> > are different. So when I manually plot a correlogram (model residuals
> > Moran's I vs distance classes), the correlogram looks different every
> > session, although I am using the same data.
>
> Unless you set the seed to the random number generator for the Monte Carlo
> procedure (set.seed()), each run will get a different stream of numbers,
> hence give possibly somewhat different results, won't it? However, the
> "statistic" component of the returned object that you are extracting is
> the same as that returned by moran.test(), and will not change for the
> same input data and spatial weights.
>
> Please take a longer coffee break (or your prefered beverage), think
> through what you are trying to do, read the relevant literature, make some
> simple examples, and you'll find that you've resolved most of your
> concerns. You are trying to do too much, too fast (deadline?), without
> grasping the essence. If you simplify, you'll see how to do things.
>
> Have you looked at correlog() in pgirmess, which effectively does what you
> "seem" to be trying to do - I pointed this out in a reply several days
> ago?
>
> Roger
>
> PS. Using moran.mc() or moran.test() on residuals is in any case
> inappropriate, you probably ought to be using lm.morantest() on the lm()
> or glm() output object (glm() not well founded in the literature).
>
> >
> > Is this normal? Or is there anything wrong with my code or neighbour weight
> > lists? When I try removing the glist, and just use (ME0_200.listw <-
> nb2listw
> > (nb0_200, style="W", zero.policy=T), I continue to get different Moran's I
> > values every time I run it.
> >
> > Here's my code:
> > ## neighbourhood bins for correlogram (0-5000km, in 200 km increments)
> > nb0_200 <- dnearneigh(as.matrix(dat$x_long, dat$y_lat), 0, 200, longlat=T)
> > nb0_200_dists <- nbdists(nb0_200, as.matrix(dat$x_long, dat$y_lat))
> > nb0_200_sims <- lapply(nb0_200_dists, function(x) (1/x))
> > ME0_200.listw <- nb2listw(nb0_200, glist=nb0_200_sims, style="W",
> zero.policy=T)
> >
> > nb200_400 <- dnearneigh(as.matrix(dat$x_long, dat$y_lat), 200, 400,
> longlat=T)
> > nb200_400_dists <- nbdists(nb200_400, as.matrix(dat$x_long, dat$y_lat))
> > nb200_400_sims <- lapply(nb200_400_dists, function(x) (1/x))
> > ME200_400.listw <- nb2listw(nb200_400, glist=nb200_400_sims, style="W",
> > zero.policy=T)
> >
> > ...
> >
> > nb4800_5000 <- dnearneigh(as.matrix(dat$x_long, dat$y_lat), 4800, 5000,
> > longlat=T)
> > nb4800_5000_dists <- nbdists(nb4800_5000, as.matrix(dat$x_long, dat$y_lat))
> > nb4800_5000_sims <- lapply(nb4800_5000_dists, function(x) (1/x))
> > ME4800_5000.listw <- nb2listw(nb4800_5000, glist=nb4800_5000_sims,
> style="W",
> > zero.policy=T)
> >
> > # Check for spatial autocorrelation using permutation tests for Moran's I
> > (Monte Carlo sim of 10000 runs)
> > n1<- moran.mc(residuals(mod1), nsim=9999, listw=ME0_200.listw,
> > na.action=na.omit, zero.policy=T)
> > n2<- moran.mc(residuals(mod1), nsim=9999, listw=ME200_400.listw,
> > na.action=na.omit, zero.policy=T)
> > ...
> >
> > n25<- moran.mc(residuals(mod1), nsim=9999, listw=ME4800_5000.listw,
> > na.action=na.omit, zero.policy=T)
> >
> > # plot correlogram of glm residuals
> > n.vec <- c(n1$statistic, n2$statistic, n3$statistic, n4$statistic, n5
> > $statistic, n6$statistic,
> > n7$statistic, n8$statistic, n9$statistic, n10$statistic, n11$statistic, n12
> > $statistic,
> > n13$statistic, n14$statistic, n15$statistic, n16$statistic, n17$statistic,
> n18
> > $statistic,
> > n19$statistic, n20$statistic, n21$statistic, n22$statistic, n23$statistic,
> n24
> > $statistic,
> > n25$statistic)
> >
> > plot(n.vec, type="b", col="red", ylim=c(-0.5,0.5))
> >
> > Many thanks,
> > Xingli
> >
> > _______________________________________________
> > R-sig-Geo mailing list
> > R-sig-Geo at stat.math.ethz.ch
> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >
>
> --
> 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
>
>