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

[Roger Bivand]

On Mon, 19 Dec 2011, Henri Vallès wrote:

> 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 choice of alpha=0.05 is very conservative, and could well be larger. 
Your results are down to your data:

library(spdep)
example(columbus)
lagcol <- SpatialFiltering(CRIME ~ 1, data=columbus, nb=col.gal.nb,
   style="W", alpha=0.05, verbose=TRUE)
lmlag <- lm(CRIME ~ fitted(lagcol), data=columbus)
lm.morantest(lmlag, nb2listw(col.gal.nb))
lagcol1 <- ME(CRIME ~ 1, data=columbus, family="gaussian",
   listw=nb2listw(col.gal.nb), alpha=0.05, stdev=TRUE, verbose=TRUE)
lmlag1 <- lm(CRIME ~ fitted(lagcol1), data=columbus)
lm.morantest(lmlag1, nb2listw(col.gal.nb))

which all show lm.morantest() at approximately the target alpha. Note that 
ME() with the default Gaussian family is the same as SpatialFiltering() 
but using a bootstrap test on the regression residuals (an approximation), 
rather than the analytical test, so SpatialFiltering() would normally be 
preferred. There is something in your data that is driving this. Please 
try to reproduce your problem with a data set shipped with spdep or 
maptools. Your null model is also asking for trouble, as most of the 
observed autocorrelation in dens is caused by missing covariates.

Hope this helps,

Roger


>
> 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
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

-- 
Roger Bivand
Department of Economics, NHH Norwegian School of Economics,
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