[R-sig-Geo] Why PCNM (MEM) approach cannot remove spatial autocorrelation in the residuals?

niv de malach nivdemalach at gmail.com
Wed Nov 15 16:31:14 CET 2017 

Now everything is clear.
Thanks a lot,
Niv
ᐧ

On Wed, Nov 15, 2017 at 1:59 PM, Roger Bivand <Roger.Bivand at nhh.no> wrote:

> On Wed, 15 Nov 2017, niv de malach wrote:
>
> Attaching the code and the dataset:
>>
>>
> Thanks for the example, very helpful:
>
>
>> # code
>> library(spdep)
>> library(adespatial)
>> mydata=read.csv("example_data.csv",header=TRUE)
>> attach(mydata)
>>
>>
> Don't use attach, in modelling use data=.
>
>
>> # creating neighbors map from the coordinates
>> xy=SpatialPoints(as.matrix(mydata[,c(5,6)]))
>> myneighbors=dnearneigh(as.matrix(mydata[,c(5,6)]), 0, 1000, row.names =
>> NULL, longlat = TRUE)# 1000 kilometer distance
>> myweight=nb2listw(myneighbors,style="W") # style W
>> # creating PCNM
>> spatial_eigenvectors=dbmem(as.matrix(mydata[,c(5,6)]),MEM.au
>> tocor="positive")
>>
>
> The body is buried here ^^^
>
> adespatial::dbmem constructs a completely different weights matrix
> internally from your coordinates. I see:
>
> attr(spatial_eigenvectors, "listw")
>>
> Characteristics of weights list object:
> Neighbour list object:
> Number of regions: 233
> Number of nonzero links: 25604
> Percentage nonzero weights: 47.16241
> Average number of links: 109.8884
>
> Weights style: B
> Weights constants summary:
>     n    nn       S0       S1       S2
> B 233 54289 25169.11 49501.47 12272171
>
> retreived by setting store.listw=TRUE - this does not respect longlat
> either, I believe, so the eigenvectors are taken from this. Your weights
> are:
>
> myweight
>>
> Characteristics of weights list object:
> Neighbour list object:
> Number of regions: 233
> Number of nonzero links: 8260
> Percentage nonzero weights: 15.21487
> Average number of links: 35.45064
>
> Weights style: W
> Weights constants summary:
>     n    nn  S0       S1       S2
> W 233 54289 233 34.27034 939.7534
>
> Using the correct lm.morantest():
>
> lm.morantest(PCNM_lm, attr(spatial_eigenvectors, "listw"),
>>
> alternative="two.sided") # significant spatial autocorelation
>
>         Global Moran I for regression residuals
>
> data:
> model: lm(formula = alpha ~ evenness + I_A + gamma + herb_cover +
> woody_cover + spatial_eigenvectors[, 1] + spatial_eigenvectors[, 2] +
> spatial_eigenvectors[, 3] + spatial_eigenvectors[, 4] +
> spatial_eigenvectors[, 5] + spatial_eigenvectors[, 6] +
> spatial_eigenvectors[, 7] + spatial_eigenvectors[, 8] +
> spatial_eigenvectors[, 9] + spatial_eigenvectors[, 10] +
> spatial_eigenvectors[, 11] + spatial_eigenvectors[, 12], data = mydata)
> weights: attr(spatial_eigenvectors, "listw")
>
> Moran I statistic standard deviate = -9.7467, p-value < 2.2e-16
> alternative hypothesis: two.sided
> sample estimates:
> Observed Moran I      Expectation         Variance
>    -2.481095e-02    -1.159734e-02     1.837913e-06
>
> lm.morantest(PCNM_lm, myweight, alternative="two.sided")
>>
>         Global Moran I for regression residuals
>
> data:
> model: lm(formula = alpha ~ evenness + I_A + gamma + herb_cover +
> woody_cover + spatial_eigenvectors[, 1] + spatial_eigenvectors[, 2] +
> spatial_eigenvectors[, 3] + spatial_eigenvectors[, 4] +
> spatial_eigenvectors[, 5] + spatial_eigenvectors[, 6] +
> spatial_eigenvectors[, 7] + spatial_eigenvectors[, 8] +
> spatial_eigenvectors[, 9] + spatial_eigenvectors[, 10] +
> spatial_eigenvectors[, 11] + spatial_eigenvectors[, 12], data = mydata)
> weights: myweight
>
> Moran I statistic standard deviate = 7.1459, p-value = 8.943e-13
> alternative hypothesis: two.sided
> sample estimates:
> Observed Moran I      Expectation         Variance
>     0.0900654150    -0.0350231024     0.0003064253
>
> The actual weights used in dbmem show significant negative residual
> autocorrelation, using your weights, you see positive, because the weights
> were not the same. Using spdep::SpatialFiltering() - see reference there -
> the eigenvectors are chosen by brute force search to remove autocorrelation
> from the residuals:
>
> SF <- SpatialFiltering(alpha ~ evenness + I_A + gamma + herb_cover +
>>
> woody_cover, data = mydata, nb=myneighbors, style="W", ExactEV=TRUE)
>
> SF_fit <- lm(alpha ~ evenness + I_A + gamma + herb_cover + woody_cover +
>>
> fitted(SF), data = mydata)
>
>> lm.morantest(SF_fit, myweight)
>>
>
>         Global Moran I for regression residuals
>
> data:
> model: lm(formula = alpha ~ evenness + I_A + gamma + herb_cover +
> woody_cover + fitted(SF), data = mydata)
> weights: myweight
>
> Moran I statistic standard deviate = 0.19667, p-value = 0.422
> alternative hypothesis: greater
> sample estimates:
> Observed Moran I      Expectation         Variance
>    -0.0333289800    -0.0370711350     0.0003620337
>
>> lm.morantest(SF_fit, attr(spatial_eigenvectors, "listw"))
>>
>
>         Global Moran I for regression residuals
>
> data:
> model: lm(formula = alpha ~ evenness + I_A + gamma + herb_cover +
> woody_cover + fitted(SF), data = mydata)
> weights: attr(spatial_eigenvectors, "listw")
>
> Moran I statistic standard deviate = 1.371, p-value = 0.08519
> alternative hypothesis: greater
> sample estimates:
> Observed Moran I      Expectation         Variance
>    -5.700770e-03    -8.090576e-03     3.038630e-06
>
> Choosing eigenvectors from the listw object returned by dbmem is rather
> time-consuming, suggesting that dbmem should pass through longlat= to be
> used internally by spdep::dnearneigh() and spdep::nbdists() used
> internally; you could then use the threshold argument. Do use
> spdep::lm.morantest() on regression residuals; note that there is no proper
> test for residual autocorrelation for fitted spatial regression models.
>
> Hope this clarifies,
>
> Roger
>
>
>> # simple regression model
>> alpha_lm=lm(alpha~evenness+I_A+gamma+herb_cover+woody_cover) #
>> summary(alpha_lm)
>> moran.test(alpha_lm$residuals,myweight) # significant spatial
>> autocorelation
>>
>>
>> # regression model with all positive eigenvectors
>> PCNM_lm=lm(alpha~evenness+I_A+gamma+herb_cover+woody_cover+
>>
>> spatial_eigenvectors[,1]+spatial_eigenvectors[,2]+spatial_
>> eigenvectors[,3]+
>>
>> spatial_eigenvectors[,4]+spatial_eigenvectors[,5]+spatial_
>> eigenvectors[,6]+
>>
>> spatial_eigenvectors[,7]+spatial_eigenvectors[,8]+spatial_
>> eigenvectors[,9]+
>>
>> spatial_eigenvectors[,10]+spatial_eigenvectors[,11]+spatial_
>> eigenvectors[,12])
>>
>> summary(PCNM_lm)
>> moran.test(PCNM_lm$residuals,myweight) # significant spatial
>> autocorelation
>>
>>
>> # SAR error regression model
>> sar_alpha=errorsarlm(alpha~evenness+I_A+gamma+herb_cover+woo
>> dy_cover,data=mydata,myweight)
>> summary(sar_alpha)
>> moran.test(sar_alpha$residuals,myweight) # no spatial autocorelation
>>
>>
>>
>>
>>
>>
>> On Wed, Nov 15, 2017 at 12:24 PM, Roger Bivand <Roger.Bivand at nhh.no>
>> wrote:
>>
>> On Wed, 15 Nov 2017, niv de malach wrote:
>>>
>>> Hi,
>>>
>>>> I am trying to include spatial eigenvectors ( in my regression using
>>>> *dbmem
>>>> *command from *adespatial* package) in order to account for spatial
>>>> correlation. The problem is that even after including all the positive
>>>> eigenvectors there is still a positive significant spatial
>>>> autocorrelation
>>>> in the residuals (based on Moran's I test). The magnitude of this
>>>> problem
>>>> is affected by the styles I use for the spatial weight (using the styles
>>>> "U","W" "B" "C" of the function *nb2list* from *spdep*) but in all
>>>> styles
>>>> Moran's I is still significantly positive.
>>>>
>>>>
>>> Maybe there is negative residual autocorrelation, and only choosing
>>> eigenvectors matching positive eigenvalues is enhancing pattern jumble
>>> (oversmoothing the model leaving jumble in the residuals)? Do you have an
>>> example you could share (link to offline downloadable code+data if need
>>> be)?
>>>
>>> Roger
>>>
>>>
>>> Interestingly this problem doesn't occur when I use a SAR models (
>>>> *errorsarlm* command from *spdep* package).
>>>>
>>>> So, does it make sense to use PCNM approach when it removes only a
>>>> portion
>>>> the spatial autocorrelation?
>>>>
>>>> Thanks
>>>> Niv
>>>>
>>>> ᐧ
>>>>
>>>>         [[alternative HTML version deleted]]
>>>>
>>>> _______________________________________________
>>>> 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, No
>>> <https://maps.google.com/?q=Economics,+No&entry=gmail&source=g>rwegian
>>> School of Economics,
>>> Helleveien 30, N-5045 Bergen, Norway.
>>> voice: +47 55 95 93 55; e-mail: Roger.Bivand at nhh.no
>>> Editor-in-Chief of The R Journal, https://journal.r-project.org/
>>> index.html
>>> http://orcid.org/0000-0003-2392-6140
>>> https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en
>>>
>>
>>
>> ᐧ
>>
>>
> --
> Roger Bivand
> Department of Economics, Norwegian School of Economics,
> Helleveien 30, N-5045 Bergen, Norway.
> voice: +47 55 95 93 55; e-mail: Roger.Bivand at nhh.no
> Editor-in-Chief of The R Journal, https://journal.r-project.org/index.html
> http://orcid.org/0000-0003-2392-6140
> https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en
>

	[[alternative HTML version deleted]]

More information about the R-sig-Geo mailing list