[R-sig-Geo] poly2nb neighbour itself should be considered a neighbour

Mon Nov 4 11:32:58 CET 2019

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 <r-sig-geo-bounces using r-project.org> on behalf of Elias T. Krainski <eliaskrainski using gmail.com>
> Sent: Sunday, November 3, 2019 23:36
> To: r-sig-geo using r-project.org <r-sig-geo using r-project.org>
> 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 <usercatch using outlook.com> 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
>>> R-sig-Geo using r-project.org
>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>        [[alternative HTML version deleted]]
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> R-sig-Geo using r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo using r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
> 	[[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo using r-project.org
> 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: Roger.Bivand using nhh.no
https://orcid.org/0000-0003-2392-6140
https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en