[R-sig-Geo] Why PCNM (MEM) approach cannot remove spatial autocorrelation in the residuals?
Roger Bivand
Roger.Bivand at nhh.no
Wed Nov 15 12:59:27 CET 2017
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.autocor="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+woody_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
More information about the R-sig-Geo
mailing list