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

[Henri Vallès]

Dear list,

Would it be appropriate to use the SpatialFiltering function of spdep
to model the spatial structure of a single variable of interest? Here,
my final objective is to carry out a variance partitioning separating
environmental "effects" from spatial "effects" and assessing the
variance in the response variable that is shared by both the space and
the environment. In this context, what I have often seen in the
literature is the use of forward selection of Moran Eigenvector Maps
to model the spatial component of the response variable, which (I
understand) is based on a different selection criterion of MEMs from
that used by SpatialFiltering.

Henri

On Tue, Dec 20, 2011 at 4:36 AM, Roger Bivand <Roger.Bivand at nhh.no> wrote:
> 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