[R-sig-Geo] removing spatial autocorrelation with ME function and subsequently testing with lm.morantest

Henri Vallès hevals at gmail.com
Mon Dec 19 15:07:43 CET 2011


Hi,

I am hoping somebody will be able to assist me.

I am using the "ME" function of the "spdep" package to select a suite
of eigenvectors that successfully removes any significant spatial
autocorrelation (at a nominal level of 0.05) from my dependent
variable (density, "den") given a spatial weighting matrix
(listwdist).

The code and output:

> ME_model=ME(den~1,listw=listwdist,verbose=TRUE,nsim=999,alpha=0.05)

eV[,1], I: 0.2267349 ZI: NA, pr(ZI): 0.001
eV[,8], I: 0.1727207 ZI: NA, pr(ZI): 0.001
eV[,3], I: 0.1490887 ZI: NA, pr(ZI): 0.001
eV[,5], I: 0.1260945 ZI: NA, pr(ZI): 0.001
eV[,2], I: 0.1042690 ZI: NA, pr(ZI): 0.001
eV[,15], I: 0.08287102 ZI: NA, pr(ZI): 0.001
eV[,10], I: 0.06740037 ZI: NA, pr(ZI): 0.001
eV[,7], I: 0.0513471 ZI: NA, pr(ZI): 0.001
eV[,21], I: 0.03391224 ZI: NA, pr(ZI): 0.018
eV[,13], I: 0.01831637 ZI: NA, pr(ZI): 0.073


However, when I use the "lm.morantest" function to double-check
whether there is any significant autocorrelation in the residuals,
using "den" as the dependent variable and the selected eigenvectors as
predictors (ME_model$vectors), the test indicates that there still is
considerable autocorrelation left in the residuals:

 model=lm(den~ME_model$vectors)

> lm.morantest(model,listwdist)

        Global Moran's I for regression residuals

data:
model: lm(formula = den ~ ME_model$vectors)
weights: listwdist

Moran I statistic standard deviate = 5.2239, p-value = 8.761e-08
alternative hypothesis: greater
sample estimates:
Observed Moran's I        Expectation           Variance
      1.831637e-02      -2.453698e-02       6.729546e-05


In contrast, I get a different result (but one that would be more
consistent with the outcome of the ME function) when I extract the
residuals from the model above and run a "moran.test" on those
residuals, although I have just found out that using the residuals
this way appears not to be valid:

>moran.test(residuals(model),listwdist)

Moran's I test under randomisation, .
data:  residuals(model)
weights: listwdist
Moran I statistic standard deviate = 1.5546, p-value = 0.06002
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance
     0.0183163706     -0.0024154589      0.0001778361


I would very much appreciate if anyone could help me understand these
discrepancies between methods or what is it that I am doing wrong?

Regards,

Henri Valles
Lecturer in Ecology
Department of Biological and Chemical Sciences
University of the West Indies, Cave Hill Campus, Barbados, BB11000
tel (Direct): (246) 417-4328
E-mail: hevals at gmail.com



More information about the R-sig-Geo mailing list