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. 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]] _______________________________________________ RsigGeo mailing list [hidden email] https://stat.ethz.ch/mailman/listinfo/rsiggeo 
Administrator

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]] > > _______________________________________________ > RsigGeo mailing list > [hidden email] > https://stat.ethz.ch/mailman/listinfo/rsiggeo Roger Bivand Department of Economics, Norwegian School of Economics, Helleveien 30, N5045 Bergen, Norway. voice: +47 55 95 93 55; email: [hidden email] EditorinChief of The R Journal, https://journal.rproject.org/index.html http://orcid.org/0000000323926140 https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en _______________________________________________ RsigGeo mailing list [hidden email] https://stat.ethz.ch/mailman/listinfo/rsiggeo
Roger Bivand
Department of Economics Norwegian School of Economics Helleveien 30 N5045 Bergen, Norway 
Attaching the code and the dataset: # code library(spdep) library(adespatial) mydata=read.csv("example_data.csv",header=TRUE) attach(mydata) # 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") # 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 <[hidden email]> wrote: On Wed, 15 Nov 2017, niv de malach wrote: ᐧ
_______________________________________________ RsigGeo mailing list [hidden email] https://stat.ethz.ch/mailman/listinfo/rsiggeo example_data.csv (29K) Download Attachment 
Administrator

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, pvalue < 2.2e16 alternative hypothesis: two.sided sample estimates: Observed Moran I Expectation Variance 2.481095e02 1.159734e02 1.837913e06 > 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, pvalue = 8.943e13 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, pvalue = 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, pvalue = 0.08519 alternative hypothesis: greater sample estimates: Observed Moran I Expectation Variance 5.700770e03 8.090576e03 3.038630e06 Choosing eigenvectors from the listw object returned by dbmem is rather timeconsuming, 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 <[hidden email]> 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]] >>> >>> _______________________________________________ >>> RsigGeo mailing list >>> [hidden email] >>> https://stat.ethz.ch/mailman/listinfo/rsiggeo >>> >> >>  >> Roger Bivand >> Department of Economics, No >> <https://maps.google.com/?q=Economics,+No&entry=gmail&source=g>rwegian >> School of Economics, >> Helleveien 30, N5045 Bergen, Norway. >> voice: +47 55 95 93 55; email: [hidden email] >> EditorinChief of The R Journal, https://journal.rproject.org/index.html >> http://orcid.org/0000000323926140 >> https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en > > > ᐧ > Roger Bivand Department of Economics, Norwegian School of Economics, Helleveien 30, N5045 Bergen, Norway. voice: +47 55 95 93 55; email: [hidden email] EditorinChief of The R Journal, https://journal.rproject.org/index.html http://orcid.org/0000000323926140 https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en _______________________________________________ RsigGeo mailing list [hidden email] https://stat.ethz.ch/mailman/listinfo/rsiggeo
Roger Bivand
Department of Economics Norwegian School of Economics Helleveien 30 N5045 Bergen, Norway 
Now everything is clear.
Thanks a lot, Niv ᐧ On Wed, Nov 15, 2017 at 1:59 PM, Roger Bivand <[hidden email]> 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, pvalue < 2.2e16 > alternative hypothesis: two.sided > sample estimates: > Observed Moran I Expectation Variance > 2.481095e02 1.159734e02 1.837913e06 > > 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, pvalue = 8.943e13 > 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, pvalue = 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, pvalue = 0.08519 > alternative hypothesis: greater > sample estimates: > Observed Moran I Expectation Variance > 5.700770e03 8.090576e03 3.038630e06 > > Choosing eigenvectors from the listw object returned by dbmem is rather > timeconsuming, 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 <[hidden email]> >> 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]] >>>> >>>> _______________________________________________ >>>> RsigGeo mailing list >>>> [hidden email] >>>> https://stat.ethz.ch/mailman/listinfo/rsiggeo >>>> >>>> >>>  >>> Roger Bivand >>> Department of Economics, No >>> <https://maps.google.com/?q=Economics,+No&entry=gmail&source=g>rwegian >>> School of Economics, >>> Helleveien 30, N5045 Bergen, Norway. >>> voice: +47 55 95 93 55; email: [hidden email] >>> EditorinChief of The R Journal, https://journal.rproject.org/ >>> index.html >>> http://orcid.org/0000000323926140 >>> https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en >>> >> >> >> ᐧ >> >> >  > Roger Bivand > Department of Economics, Norwegian School of Economics, > Helleveien 30, N5045 Bergen, Norway. > voice: +47 55 95 93 55; email: [hidden email] > EditorinChief of The R Journal, https://journal.rproject.org/index.html > http://orcid.org/0000000323926140 > https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en > [[alternative HTML version deleted]] _______________________________________________ RsigGeo mailing list [hidden email] https://stat.ethz.ch/mailman/listinfo/rsiggeo 
Free forum by Nabble  Edit this page 