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 <-
(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",
nb200_400 <- dnearneigh(as.matrix(dat$x_long, dat$y_lat), 200, 400,
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,
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,
$statistic,
n19$statistic, n20$statistic, n21$statistic, n22$statistic, n23$statistic,
$statistic,
n25$statistic)
plot(n.vec, type="b", col="red", ylim=c(-0.5,0.5))
Many thanks,
Xingli