poly2nb neighbour itself should be considered a neighbour

classic Classic list List threaded Threaded
8 messages Options
Reply | Threaded
Open this post in threaded view
|

poly2nb neighbour itself should be considered a neighbour

Robert R
Dear All,

I would like to know if the function "poly2nb" ("spdep" package.) let me create a neighborhood of itself, i.e., not only its queen neighbors (queen=TRUE), but a neighbour itself should also be considered a neighbour.

I am looking to create a queen weight neighborhood matrix afterwards using "nb2listw".

Any help would help me a lot.

Many thanks

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: poly2nb neighbour itself should be considered a neighbour

Dexter Locke
Dear Robert,

It sounds like what you are looking for is typically called a second order neighbor. Higher order neighbors can also included in a weights matrix such as your neighbors’, neighbors’, neighbor which is a third-order neighbor. I think you are seeking second-order neighbors.

See the spdep vignettes, and the section 5 Higher-Order Neighbors subsection here: https://cran.r-project.org/web/packages/spdep/vignettes/nb.pdf in particular. The spdep::nblag might be what you need, but without additional information it is hard to know.

Good luck,
Dexter
dexterlocke.com



> On Nov 3, 2019, at 2:12 PM, Robert R <[hidden email]> wrote:
>
> Dear All,
>
> I would like to know if the function "poly2nb" ("spdep" package.) let me create a neighborhood of itself, i.e., not only its queen neighbors (queen=TRUE), but a neighbour itself should also be considered a neighbour.
>
> I am looking to create a queen weight neighborhood matrix afterwards using "nb2listw".
>
> Any help would help me a lot.
>
> Many thanks
>
>    [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: poly2nb neighbour itself should be considered a neighbour

Elias T Krainski
Hello,

I found this matter easier when working with matrix representations. Set

  A^0 = I (identity matrix),

  A^1 = A, where A_{ij} = 1 if j is neighbor to j and zero otherwise
(this consider A_{ii} = 0)

  A^2 = A'A

  A^3 = A'A^2 and so on

The A^k_{ij} entry gives _how many steps of length k there is between i
and j_. To me, this definition makes this matter clear. See an example
considering a 10x10 regular grid:

nb <- grid2nb(d=c(10,10), queen = FALSE)
nn <- card(nb)
A1 <- sparseMatrix(
  i = rep(1:length(nn), nn),
  j = unlist(nb[nn>0]), x=1)
A2 <- crossprod(A1)
image(A1)
image(A2)

best regards,

Elias

On 03/11/2019 17:46, Dexter Locke wrote:

> Dear Robert,
>
> It sounds like what you are looking for is typically called a second order neighbor. Higher order neighbors can also included in a weights matrix such as your neighbors’, neighbors’, neighbor which is a third-order neighbor. I think you are seeking second-order neighbors.
>
> See the spdep vignettes, and the section 5 Higher-Order Neighbors subsection here: https://cran.r-project.org/web/packages/spdep/vignettes/nb.pdf in particular. The spdep::nblag might be what you need, but without additional information it is hard to know.
>
> Good luck,
> Dexter
> dexterlocke.com
>
>
>
>> On Nov 3, 2019, at 2:12 PM, Robert R <[hidden email]> wrote:
>>
>> Dear All,
>>
>> I would like to know if the function "poly2nb" ("spdep" package.) let me create a neighborhood of itself, i.e., not only its queen neighbors (queen=TRUE), but a neighbour itself should also be considered a neighbour.
>>
>> I am looking to create a queen weight neighborhood matrix afterwards using "nb2listw".
>>
>> Any help would help me a lot.
>>
>> Many thanks
>>
>>     [[alternative HTML version deleted]]
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> [hidden email]
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: poly2nb neighbour itself should be considered a neighbour

Robert R
Dear Dexter, Dear Elias,

Many thanks for your response.

I am looking to create a spatial weight matrix (queen contiguity neighbors), with the only difference that a neighbour itself should also be considered a neighbour.

Below I am sending you the code. Basically I create a polygon for NYC and give that polygon I want to create the spatial weight matrix.

@Elias: I used your answer to generate two matrices A1 and A2. But how to turn them into an nb object so I can use the function nb2listw (with the arguments style="W", zero.policy=TRUE).

Thank you and best regards,
Robert


#####

####################
# ---- packages ----​
####################​

packages_install <- function(packages){​
 new.packages <- packages[!(packages %in% installed.packages()[, "Package"])]​
 if (length(new.packages)) ​
 install.packages(new.packages, dependencies = TRUE)​
 sapply(packages, require, character.only = TRUE)​
}​

packages_required <- c("data.table", "dplyr", "sf", "spdep", "Matrix")​
packages_install(packages_required)​

# Working directory​
setwd("C:/Users/User/Documents/Code")​


#####################​
# ---- zips_nyc ----​
#####################​

zips_nyc_bronx <- c("10451", "10452", "10453", "10454", "10455", "10456", "10457", "10458", "10459", "10460", "10461", "10462", "10463", "10464", "10465", "10466", "10467", "10468", "10469", "10470", "10471", "10472", "10473", "10474", "10475")​
zips_nyc_brooklyn <- c("11201", "11203", "11204", "11205", "11206", "11207", "11208", "11209", "11210", "11211", "11212", "11213", "11214", "11215", "11216", "11217", "11218", "11219", "11220", "11221", "11222", "11223", "11224", "11225", "11226", "11228", "11229", "11230", "11231", "11232", "11233", "11234", "11235", "11236", "11237", "11238", "11239")​
zips_nyc_manhattan <- c("10001", "10002", "10003", "10004", "10005", "10006", "10007", "10009", "10010", "10011", "10012", "10013", "10014", "10016", "10017", "10018", "10019", "10020", "10021", "10022", "10023", "10024", "10025", "10026", "10027", "10028", "10029", "10030", "10031", "10032", "10033", "10034", "10035", "10036", "10037", "10038", "10039", "10040", "10044", "10065", "10075", "10128", "10280")​
zips_nyc_queens <- c("11004", "11005", "11101", "11102", "11103", "11104", "11105", "11106", "11354", "11355", "11356", "11357", "11358", "11359", "11360", "11361", "11362", "11363", "11364", "11365", "11366", "11367", "11368", "11369", "11370", "11372", "11373", "11374", "11375", "11377", "11378", "11379", "11385", "11411", "11412", "11413", "11414", "11415", "11416", "11417", "11418", "11419", "11420", "11421", "11422", "11423", "11426", "11427", "11428", "11429", "11432", "11433", "11434", "11435", "11436", "11691", "11692", "11693", "11694", "11695", "11697")​
zips_nyc_staten_island <- c("10301", "10302", "10303", "10304", "10305", "10306", "10307", "10308", "10309", "10310", "10312", "10314")​
zips_nyc <- sort(c(zips_nyc_bronx, zips_nyc_brooklyn, zips_nyc_manhattan, zips_nyc_queens, zips_nyc_staten_island))​


#####################​
# ---- shapefile ----​
#####################​

## shapefile_us​

# Shapefile zips import and Coordinate Reference System (CRS) transformation​
# Download: https://www2.census.gov/geo/tiger/GENZ2018/shp/cb_2018_us_zcta510_500k.zip​
shapefile_us <- sf::st_read(dsn = "Shapefile", layer = "cb_2018_us_zcta510_500k")​

# Columns removal​
shapefile_us <- shapefile_us %>% select(-c(AFFGEOID10, GEOID10, ALAND10, AWATER10))​

# Column rename: ZCTA5CE10​
setnames(shapefile_us, old=c("ZCTA5CE10"), new=c("zipcode"))​

# Column class change: zipcode​
shapefile_us$zipcode <- as.character(shapefile_us$zipcode)​


## polygon_nyc​
polygon_nyc <- shapefile_us %>% filter(zipcode %in% zips_nyc)​

# Variable creation: list of neighbors for each polygon (queen contiguity neighbors)​
nb <- poly2nb(polygon_nyc, queen=FALSE)​
nn <- card(nb)​
A1 <- sparseMatrix(​
  i = rep(1:length(nn), nn),​
  j = unlist(nb[nn>0]), x=1)​
A2 <- crossprod(A1)​
image(A1)​
image(A2)​

# next, supplement the neighbor list with spatial weights: "W" row-standardize weights​
W_Matrix <- nb2listw(neighbours = A1, style="W", zero.policy=TRUE)

#####
________________________________
From: R-sig-Geo <[hidden email]> on behalf of Elias T. Krainski <[hidden email]>
Sent: Sunday, November 3, 2019 23:36
To: [hidden email] <[hidden email]>
Subject: Re: [R-sig-Geo] poly2nb neighbour itself should be considered a neighbour

Hello,

I found this matter easier when working with matrix representations. Set

  A^0 = I (identity matrix),

  A^1 = A, where A_{ij} = 1 if j is neighbor to j and zero otherwise
(this consider A_{ii} = 0)

  A^2 = A'A

  A^3 = A'A^2 and so on

The A^k_{ij} entry gives _how many steps of length k there is between i
and j_. To me, this definition makes this matter clear. See an example
considering a 10x10 regular grid:

nb <- grid2nb(d=c(10,10), queen = FALSE)
nn <- card(nb)
A1 <- sparseMatrix(
  i = rep(1:length(nn), nn),
  j = unlist(nb[nn>0]), x=1)
A2 <- crossprod(A1)
image(A1)
image(A2)

best regards,

Elias

On 03/11/2019 17:46, Dexter Locke wrote:

> Dear Robert,
>
> It sounds like what you are looking for is typically called a second order neighbor. Higher order neighbors can also included in a weights matrix such as your neighbors’, neighbors’, neighbor which is a third-order neighbor. I think you are seeking second-order neighbors.
>
> See the spdep vignettes, and the section 5 Higher-Order Neighbors subsection here: https://cran.r-project.org/web/packages/spdep/vignettes/nb.pdf in particular. The spdep::nblag might be what you need, but without additional information it is hard to know.
>
> Good luck,
> Dexter
> dexterlocke.com
>
>
>
>> On Nov 3, 2019, at 2:12 PM, Robert R <[hidden email]> wrote:
>>
>> Dear All,
>>
>> I would like to know if the function "poly2nb" ("spdep" package.) let me create a neighborhood of itself, i.e., not only its queen neighbors (queen=TRUE), but a neighbour itself should also be considered a neighbour.
>>
>> I am looking to create a queen weight neighborhood matrix afterwards using "nb2listw".
>>
>> Any help would help me a lot.
>>
>> Many thanks
>>
>>     [[alternative HTML version deleted]]
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> [hidden email]
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>        [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: poly2nb neighbour itself should be considered a neighbour

Barry Rowlingson-3
In reply to this post by Robert R
Can you manipulate the adjacency list structure to add `i` to each list
element vector?

eg using sample data from spdep:

make a neighbour structure:

 > colnn = poly2nb(columbus)

this is a list - so for example polygon 4 is next to:

 > colnn[[4]]
 [1] 2 3 5 8

2, 3, 5, and 8. It seems you want to include `4` in that vector. So run
this loop:

 > for(i in 1:length(colnn)){colnn[[i]]=as.integer(c(i,colnn[[i]]))}

 which produces a nb structure:

 > colnn
Neighbour list object:
Number of regions: 49
Number of nonzero links: 285
Percentage nonzero weights: 11.87005
Average number of links: 5.816327

compared with the original without self-links:

> poly2nb(columbus)
Neighbour list object:
Number of regions: 49
Number of nonzero links: 236
Percentage nonzero weights: 9.829238
Average number of links: 4.816327

If you want to go on to make a weights list object:

> lw = nb2listw(colnn)

the neighbours are preserved:

> str(lw)
List of 3
 $ style     : chr "W"
 $ neighbours:List of 49
  ..$ : int [1:3] 1 2 3
  ..$ : int [1:4] 2 1 3 4
  ..$ : int [1:5] 3 1 2 4 5
  ..$ : int [1:5] 4 2 3 5 8
...
$ weights   :List of 49
  ..$ : num [1:3] 0.333 0.333 0.333
  ..$ : num [1:4] 0.25 0.25 0.25 0.25
  ..$ : num [1:5] 0.2 0.2 0.2 0.2 0.2
  ..$ : num [1:5] 0.2 0.2 0.2 0.2 0.2

So I think that's what you want...

Note carefully the use of `as.integer` here:

  for(i in 1:length(colnn)){colnn[[i]]=as.integer(c(i,colnn[[i]]))}

because there's code in `spdep` that expects these things to be stored as
integer and passes them to C code. Get this wrong when manipulating

 > xx = poly2nb(columbus)
 > xx
Neighbour list object:
Number of regions: 49
Number of nonzero links: 236
Percentage nonzero weights: 9.829238
Average number of links: 4.816327
 > xx[[1]] = c(1, xx[[1]])
 > xx
 Error in card(nb) :
  INTEGER() can only be applied to a 'integer', not a 'double'

so instead:

> xx = poly2nb(columbus)
> xx[[1]] = as.integer(c(1, xx[[1]]))
> xx
Neighbour list object:
Number of regions: 49
Number of nonzero links: 237
Percentage nonzero weights: 9.870887
Average number of links: 4.836735

works - I don't think its strictly necessary in the `for` loop because `i`
is an integer but belt and braces....

Barry



On Sun, Nov 3, 2019 at 7:12 PM Robert R <[hidden email]> wrote:

> Dear All,
>
> I would like to know if the function "poly2nb" ("spdep" package.) let me
> create a neighborhood of itself, i.e., not only its queen neighbors
> (queen=TRUE), but a neighbour itself should also be considered a neighbour.
>
> I am looking to create a queen weight neighborhood matrix afterwards using
> "nb2listw".
>
> Any help would help me a lot.
>
> Many thanks
>
>         [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: poly2nb neighbour itself should be considered a neighbour

Roger Bivand
Administrator
In reply to this post by Robert R
On Mon, 4 Nov 2019, Robert R wrote:

> Dear Dexter, Dear Elias,
>
> Many thanks for your response.
>
> I am looking to create a spatial weight matrix (queen contiguity
> neighbors), with the only difference that a neighbour itself should also
> be considered a neighbour.

I think you are looking for spdep::include.self(), originally written to
include self-neighbours in nb objects for the local Getis-Ord G_i^*
measure. Powering the matrix will also give self-neighbours as a
by-product for observations with neighbours, but not only ones on the
principal diagonal.

Roger

>
> Below I am sending you the code. Basically I create a polygon for NYC
> and give that polygon I want to create the spatial weight matrix.
>
> @Elias: I used your answer to generate two matrices A1 and A2. But how to turn them into an nb object so I can use the function nb2listw (with the arguments style="W", zero.policy=TRUE).
>
> Thank you and best regards,
> Robert
>
>
> #####
>
> ####################
> # ---- packages ----​
> ####################​
> ​
> packages_install <- function(packages){​
> new.packages <- packages[!(packages %in% installed.packages()[, "Package"])]​
> if (length(new.packages)) ​
> install.packages(new.packages, dependencies = TRUE)​
> sapply(packages, require, character.only = TRUE)​
> }​
> ​
> packages_required <- c("data.table", "dplyr", "sf", "spdep", "Matrix")​
> packages_install(packages_required)​
> ​
> # Working directory​
> setwd("C:/Users/User/Documents/Code")​
> ​
> ​
> #####################​
> # ---- zips_nyc ----​
> #####################​
> ​
> zips_nyc_bronx <- c("10451", "10452", "10453", "10454", "10455", "10456", "10457", "10458", "10459", "10460", "10461", "10462", "10463", "10464", "10465", "10466", "10467", "10468", "10469", "10470", "10471", "10472", "10473", "10474", "10475")​
> zips_nyc_brooklyn <- c("11201", "11203", "11204", "11205", "11206", "11207", "11208", "11209", "11210", "11211", "11212", "11213", "11214", "11215", "11216", "11217", "11218", "11219", "11220", "11221", "11222", "11223", "11224", "11225", "11226", "11228", "11229", "11230", "11231", "11232", "11233", "11234", "11235", "11236", "11237", "11238", "11239")​
> zips_nyc_manhattan <- c("10001", "10002", "10003", "10004", "10005", "10006", "10007", "10009", "10010", "10011", "10012", "10013", "10014", "10016", "10017", "10018", "10019", "10020", "10021", "10022", "10023", "10024", "10025", "10026", "10027", "10028", "10029", "10030", "10031", "10032", "10033", "10034", "10035", "10036", "10037", "10038", "10039", "10040", "10044", "10065", "10075", "10128", "10280")​
> zips_nyc_queens <- c("11004", "11005", "11101", "11102", "11103", "11104", "11105", "11106", "11354", "11355", "11356", "11357", "11358", "11359", "11360", "11361", "11362", "11363", "11364", "11365", "11366", "11367", "11368", "11369", "11370", "11372", "11373", "11374", "11375", "11377", "11378", "11379", "11385", "11411", "11412", "11413", "11414", "11415", "11416", "11417", "11418", "11419", "11420", "11421", "11422", "11423", "11426", "11427", "11428", "11429", "11432", "11433", "11434", "11435", "11436", "11691", "11692", "11693", "11694", "11695", "11697")​
> zips_nyc_staten_island <- c("10301", "10302", "10303", "10304", "10305", "10306", "10307", "10308", "10309", "10310", "10312", "10314")​
> zips_nyc <- sort(c(zips_nyc_bronx, zips_nyc_brooklyn, zips_nyc_manhattan, zips_nyc_queens, zips_nyc_staten_island))​
> ​
> ​
> #####################​
> # ---- shapefile ----​
> #####################​
> ​
> ## shapefile_us​
> ​
> # Shapefile zips import and Coordinate Reference System (CRS) transformation​
> # Download: https://www2.census.gov/geo/tiger/GENZ2018/shp/cb_2018_us_zcta510_500k.zip​
> shapefile_us <- sf::st_read(dsn = "Shapefile", layer = "cb_2018_us_zcta510_500k")​
> ​
> # Columns removal​
> shapefile_us <- shapefile_us %>% select(-c(AFFGEOID10, GEOID10, ALAND10, AWATER10))​
> ​
> # Column rename: ZCTA5CE10​
> setnames(shapefile_us, old=c("ZCTA5CE10"), new=c("zipcode"))​
> ​
> # Column class change: zipcode​
> shapefile_us$zipcode <- as.character(shapefile_us$zipcode)​
> ​
> ​
> ## polygon_nyc​
> polygon_nyc <- shapefile_us %>% filter(zipcode %in% zips_nyc)​
> ​
> # Variable creation: list of neighbors for each polygon (queen contiguity neighbors)​
> nb <- poly2nb(polygon_nyc, queen=FALSE)​
> nn <- card(nb)​
> A1 <- sparseMatrix(​
>  i = rep(1:length(nn), nn),​
>  j = unlist(nb[nn>0]), x=1)​
> A2 <- crossprod(A1)​
> image(A1)​
> image(A2)​
> ​
> # next, supplement the neighbor list with spatial weights: "W" row-standardize weights​
> W_Matrix <- nb2listw(neighbours = A1, style="W", zero.policy=TRUE)
>
> #####
> ________________________________
> From: R-sig-Geo <[hidden email]> on behalf of Elias T. Krainski <[hidden email]>
> Sent: Sunday, November 3, 2019 23:36
> To: [hidden email] <[hidden email]>
> Subject: Re: [R-sig-Geo] poly2nb neighbour itself should be considered a neighbour
>
> Hello,
>
> I found this matter easier when working with matrix representations. Set
>
>  A^0 = I (identity matrix),
>
>  A^1 = A, where A_{ij} = 1 if j is neighbor to j and zero otherwise
> (this consider A_{ii} = 0)
>
>  A^2 = A'A
>
>  A^3 = A'A^2 and so on
>
> The A^k_{ij} entry gives _how many steps of length k there is between i
> and j_. To me, this definition makes this matter clear. See an example
> considering a 10x10 regular grid:
>
> nb <- grid2nb(d=c(10,10), queen = FALSE)
> nn <- card(nb)
> A1 <- sparseMatrix(
>  i = rep(1:length(nn), nn),
>  j = unlist(nb[nn>0]), x=1)
> A2 <- crossprod(A1)
> image(A1)
> image(A2)
>
> best regards,
>
> Elias
>
> On 03/11/2019 17:46, Dexter Locke wrote:
>> Dear Robert,
>>
>> It sounds like what you are looking for is typically called a second order neighbor. Higher order neighbors can also included in a weights matrix such as your neighbors’, neighbors’, neighbor which is a third-order neighbor. I think you are seeking second-order neighbors.
>>
>> See the spdep vignettes, and the section 5 Higher-Order Neighbors subsection here: https://cran.r-project.org/web/packages/spdep/vignettes/nb.pdf in particular. The spdep::nblag might be what you need, but without additional information it is hard to know.
>>
>> Good luck,
>> Dexter
>> dexterlocke.com
>>
>>
>>
>>> On Nov 3, 2019, at 2:12 PM, Robert R <[hidden email]> wrote:
>>>
>>> Dear All,
>>>
>>> I would like to know if the function "poly2nb" ("spdep" package.) let me create a neighborhood of itself, i.e., not only its queen neighbors (queen=TRUE), but a neighbour itself should also be considered a neighbour.
>>>
>>> I am looking to create a queen weight neighborhood matrix afterwards using "nb2listw".
>>>
>>> Any help would help me a lot.
>>>
>>> Many thanks
>>>
>>>     [[alternative HTML version deleted]]
>>>
>>> _______________________________________________
>>> R-sig-Geo mailing list
>>> [hidden email]
>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>        [[alternative HTML version deleted]]
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> [hidden email]
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
--
Roger Bivand
Department of Economics, Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; e-mail: [hidden email]
https://orcid.org/0000-0003-2392-6140
https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en
_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Roger Bivand
Department of Economics
Norwegian School of Economics
Helleveien 30
N-5045 Bergen, Norway
Reply | Threaded
Open this post in threaded view
|

Re: poly2nb neighbour itself should be considered a neighbour

Roger Bivand
Administrator
In reply to this post by Barry Rowlingson-3
On Mon, 4 Nov 2019, Barry Rowlingson wrote:

> Can you manipulate the adjacency list structure to add `i` to each list
> element vector?
>
> eg using sample data from spdep:
>
> make a neighbour structure:
>
> > colnn = poly2nb(columbus)
>
> this is a list - so for example polygon 4 is next to:
>
> > colnn[[4]]
> [1] 2 3 5 8
>
> 2, 3, 5, and 8. It seems you want to include `4` in that vector. So run
> this loop:
>
> > for(i in 1:length(colnn)){colnn[[i]]=as.integer(c(i,colnn[[i]]))}
>
> which produces a nb structure:
>
> > colnn
> Neighbour list object:
> Number of regions: 49
> Number of nonzero links: 285
> Percentage nonzero weights: 11.87005
> Average number of links: 5.816327
>
> compared with the original without self-links:
>
>> poly2nb(columbus)
> Neighbour list object:
> Number of regions: 49
> Number of nonzero links: 236
> Percentage nonzero weights: 9.829238
> Average number of links: 4.816327
>

Barry, thanks: the include.self() function:

include.self <- function(nb) {
  if (!is.null(attributes(nb)$self.included) &&
  (as.logical(attributes(nb)$self.included)))
  stop("Self already included")
  n <- length(nb)
  nc <- card(nb)
  for (i in 1:n) {
  if (nc[i] > 0) {
  nb[[i]] <- sort(c(i, nb[[i]]))
  } else {
  nb[[i]] <- i
  }
  }

  attr(nb, "self.included") <- TRUE
  nb
}

does this with sorting, but maybe needs class inheritance checking - it
was probably written last century.

Roger

> If you want to go on to make a weights list object:
>
>> lw = nb2listw(colnn)
>
> the neighbours are preserved:
>
>> str(lw)
> List of 3
> $ style     : chr "W"
> $ neighbours:List of 49
>  ..$ : int [1:3] 1 2 3
>  ..$ : int [1:4] 2 1 3 4
>  ..$ : int [1:5] 3 1 2 4 5
>  ..$ : int [1:5] 4 2 3 5 8
> ...
> $ weights   :List of 49
>  ..$ : num [1:3] 0.333 0.333 0.333
>  ..$ : num [1:4] 0.25 0.25 0.25 0.25
>  ..$ : num [1:5] 0.2 0.2 0.2 0.2 0.2
>  ..$ : num [1:5] 0.2 0.2 0.2 0.2 0.2
>
> So I think that's what you want...
>
> Note carefully the use of `as.integer` here:
>
>  for(i in 1:length(colnn)){colnn[[i]]=as.integer(c(i,colnn[[i]]))}
>
> because there's code in `spdep` that expects these things to be stored as
> integer and passes them to C code. Get this wrong when manipulating
>
> > xx = poly2nb(columbus)
> > xx
> Neighbour list object:
> Number of regions: 49
> Number of nonzero links: 236
> Percentage nonzero weights: 9.829238
> Average number of links: 4.816327
> > xx[[1]] = c(1, xx[[1]])
> > xx
> Error in card(nb) :
>  INTEGER() can only be applied to a 'integer', not a 'double'
>
> so instead:
>
>> xx = poly2nb(columbus)
>> xx[[1]] = as.integer(c(1, xx[[1]]))
>> xx
> Neighbour list object:
> Number of regions: 49
> Number of nonzero links: 237
> Percentage nonzero weights: 9.870887
> Average number of links: 4.836735
>
> works - I don't think its strictly necessary in the `for` loop because `i`
> is an integer but belt and braces....
>
> Barry
>
>
>
> On Sun, Nov 3, 2019 at 7:12 PM Robert R <[hidden email]> wrote:
>
>> Dear All,
>>
>> I would like to know if the function "poly2nb" ("spdep" package.) let me
>> create a neighborhood of itself, i.e., not only its queen neighbors
>> (queen=TRUE), but a neighbour itself should also be considered a neighbour.
>>
>> I am looking to create a queen weight neighborhood matrix afterwards using
>> "nb2listw".
>>
>> Any help would help me a lot.
>>
>> Many thanks
>>
>>         [[alternative HTML version deleted]]
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> [hidden email]
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

--
Roger Bivand
Department of Economics, Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; e-mail: [hidden email]
https://orcid.org/0000-0003-2392-6140
https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Roger Bivand
Department of Economics
Norwegian School of Economics
Helleveien 30
N-5045 Bergen, Norway
Reply | Threaded
Open this post in threaded view
|

Re: poly2nb neighbour itself should be considered a neighbour

Robert R
Dear All,

Many thanks for the proposed solutions, specially to Roger and Barry.

This was exactly what I was looking for.

I tested with a really simple polygon (see below), and it worked.

##
columbus <- sf::read_sf(system.file("etc/shapes/columbus.shp", package="spdep"))
columbus <- columbus %>% slice(1:2)​
# plot(columbus)​

colnn <- poly2nb(columbus)​
colnn[1]​
colnn[2]​

# Solution 1​
for(i in 1:length(colnn)){colnn[[i]]=as.integer(c(i,colnn[[i]]))}​
colnn[1]​
colnn[2]​

# Solution 2​
colnn <- include.self(colnn)​
colnn[1]​
colnn[2]​

lw <- nb2listw(neighbours = colnn, style="W", zero.policy=TRUE)

##

Best regards,
Robert
________________________________
From: Roger Bivand <[hidden email]>
Sent: Monday, November 4, 2019 11:37
To: Barry Rowlingson <[hidden email]>
Cc: Robert R <[hidden email]>; [hidden email] <[hidden email]>
Subject: Re: [R-sig-Geo] poly2nb neighbour itself should be considered a neighbour

On Mon, 4 Nov 2019, Barry Rowlingson wrote:

> Can you manipulate the adjacency list structure to add `i` to each list
> element vector?
>
> eg using sample data from spdep:
>
> make a neighbour structure:
>
> > colnn = poly2nb(columbus)
>
> this is a list - so for example polygon 4 is next to:
>
> > colnn[[4]]
> [1] 2 3 5 8
>
> 2, 3, 5, and 8. It seems you want to include `4` in that vector. So run
> this loop:
>
> > for(i in 1:length(colnn)){colnn[[i]]=as.integer(c(i,colnn[[i]]))}
>
> which produces a nb structure:
>
> > colnn
> Neighbour list object:
> Number of regions: 49
> Number of nonzero links: 285
> Percentage nonzero weights: 11.87005
> Average number of links: 5.816327
>
> compared with the original without self-links:
>
>> poly2nb(columbus)
> Neighbour list object:
> Number of regions: 49
> Number of nonzero links: 236
> Percentage nonzero weights: 9.829238
> Average number of links: 4.816327
>

Barry, thanks: the include.self() function:

include.self <- function(nb) {
         if (!is.null(attributes(nb)$self.included) &&
                 (as.logical(attributes(nb)$self.included)))
                 stop("Self already included")
         n <- length(nb)
         nc <- card(nb)
         for (i in 1:n) {
                 if (nc[i] > 0) {
                         nb[[i]] <- sort(c(i, nb[[i]]))
                 } else {
                         nb[[i]] <- i
                 }
         }

         attr(nb, "self.included") <- TRUE
         nb
}

does this with sorting, but maybe needs class inheritance checking - it
was probably written last century.

Roger

> If you want to go on to make a weights list object:
>
>> lw = nb2listw(colnn)
>
> the neighbours are preserved:
>
>> str(lw)
> List of 3
> $ style     : chr "W"
> $ neighbours:List of 49
>  ..$ : int [1:3] 1 2 3
>  ..$ : int [1:4] 2 1 3 4
>  ..$ : int [1:5] 3 1 2 4 5
>  ..$ : int [1:5] 4 2 3 5 8
> ...
> $ weights   :List of 49
>  ..$ : num [1:3] 0.333 0.333 0.333
>  ..$ : num [1:4] 0.25 0.25 0.25 0.25
>  ..$ : num [1:5] 0.2 0.2 0.2 0.2 0.2
>  ..$ : num [1:5] 0.2 0.2 0.2 0.2 0.2
>
> So I think that's what you want...
>
> Note carefully the use of `as.integer` here:
>
>  for(i in 1:length(colnn)){colnn[[i]]=as.integer(c(i,colnn[[i]]))}
>
> because there's code in `spdep` that expects these things to be stored as
> integer and passes them to C code. Get this wrong when manipulating
>
> > xx = poly2nb(columbus)
> > xx
> Neighbour list object:
> Number of regions: 49
> Number of nonzero links: 236
> Percentage nonzero weights: 9.829238
> Average number of links: 4.816327
> > xx[[1]] = c(1, xx[[1]])
> > xx
> Error in card(nb) :
>  INTEGER() can only be applied to a 'integer', not a 'double'
>
> so instead:
>
>> xx = poly2nb(columbus)
>> xx[[1]] = as.integer(c(1, xx[[1]]))
>> xx
> Neighbour list object:
> Number of regions: 49
> Number of nonzero links: 237
> Percentage nonzero weights: 9.870887
> Average number of links: 4.836735
>
> works - I don't think its strictly necessary in the `for` loop because `i`
> is an integer but belt and braces....
>
> Barry
>
>
>
> On Sun, Nov 3, 2019 at 7:12 PM Robert R <[hidden email]> wrote:
>
>> Dear All,
>>
>> I would like to know if the function "poly2nb" ("spdep" package.) let me
>> create a neighborhood of itself, i.e., not only its queen neighbors
>> (queen=TRUE), but a neighbour itself should also be considered a neighbour.
>>
>> I am looking to create a queen weight neighborhood matrix afterwards using
>> "nb2listw".
>>
>> Any help would help me a lot.
>>
>> Many thanks
>>
>>         [[alternative HTML version deleted]]
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> [hidden email]
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>
>        [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

--
Roger Bivand
Department of Economics, Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; e-mail: [hidden email]
https://orcid.org/0000-0003-2392-6140
https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo