Skip to content

Different Moran's I for same model residuals in different runs using moran.mc{spdep}

2 messages · Xingli Giam, Roger Bivand

Original
Xingli Giam
# 
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.

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
Roger Bivand
#
On Sat, 31 Jan 2009, Xingli Giam wrote:

            

      
      
      
    

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