poly2nb provide 0 neighbor in US county shapefile

[Vickyxubba]

Hello everyone,

I am trying to define number of neighbor for 3108 US counties using poly2nb, however, I always have three counties have "0" neighbors... I have try to change the value of snap in poly2nb function, but without luck. Could anyone help me with this issue? Thanks in advance!

Here is the link that I download the data, the US 3143 county shapefile.

https://hub.arcgis.com/datasets/48f9af87daa241c4b267c5931ad3b226/explore

bCRS<-CRS("+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=37.5 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs")

ustry<-readShapePoly("USA_Counties.shp", proj4string=bCRS)

us.3108=ustry[ustry$STATE_FIPS != '02' & ustry$STATE_FIPS != '15', ] ---no alaska and hawaii

us.poly<-poly2nb(us.3108, snap=1000, queen=FALSE) have tried different values of snap, but no luck :/

us.nb<-nb2WB(us.poly)

us.nb$num

the num gave me three counties with "0" neighbors; however, I suppose every county should have at least one neighbors.

[rsbivand]

Please add the output of sessionInfo(), as package versions may matter.

[Vickyxubba]

Here is the output. Thank you so much.

[rsbivand]

No, the complete output, and copied as plain text, not a screen dump. Your extract does not show anything beyond loaded packages, excluding spdep itself, nor your R version and platform.

[Vickyxubba]

Sorry, I misunderstood...Do you mean something like this?

R version 4.1.2 (2021-11-01) Platform: x86_64-apple-darwin17.0 (64-bit) Running under: macOS Big Sur 11.7.2

Matrix products: default LAPACK: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRlapack.dylib

locale: [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages: character(0)

other attached packages: [1] spdep_1.1-12

[rsbivand]

Yes, of course. Next, upgrade R to current, and update all the packages. Nobody can reproduce issues found in outdated installations, developers only run updated systems.

[Vickyxubba]

I upgraded R to current and updated all package. However, I still have those three counties has 0 neighbors...

R version 4.3.1 (2023-06-16) Platform: x86_64-apple-darwin20 (64-bit) Running under: macOS Big Sur 11.7.2

Matrix products: default BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib LAPACK: /Library/Frameworks/R.framework/Versions/4.3-x86_64/Resources/lib/libRlapack.dylib; LAPACK version 3.11.0

locale: [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

time zone: America/New_York tzcode source: internal

attached base packages: character(0)

other attached packages: [1] spdep_1.2-8

[rsbivand]

1. Note that the file to be downloaded is deprecated by ESRI, so best to use an equivalent US Census file, for example through the tidycensus package.

2. The ESRI shapefile is in geographical coordinates, not in the CRS you assign. If you wanted to transform to that CRS, you must do that as a separate step.

3: my analysis:

library(sf)
library(spdep)
ustry <- st_read("USA_Counties.shp")
file.show("USA_Counties.prj")
st_crs("ESRI:102003")
st_crs("ESRI:102003")$proj4string
ustry1 <- st_transform(ustry, "ESRI:102003")
us.3108 <- ustry1[ustry1$STATE_FIPS != '02' & ustry1$STATE_FIPS != '15', ]
us.3108
us.poly_q <- poly2nb(us.3108, queen=FALSE)
us.poly <- poly2nb(us.3108)
us.poly
us.poly_q
us.3108[card(us.poly_q) == 0L, c("NAME", "FIPS"), drop=TRUE]

For these three counties: https://en.wikipedia.org/wiki/Dukes_County,_Massachusetts, https://en.wikipedia.org/wiki/Nantucket, and https://en.wikipedia.org/wiki/San_Juan_County,_Washington, the chosen contiguity criteria work as specified - the observations are real islands.

Since you want to use WinBUGS, and CAR models in general do not accept observations with no neighbours, you could choose simply to omit those observations from the whole analysis, leaving 3105 counties. You may also add in the nearest coastal county:

us_k <- knn2nb(knearneigh(st_centroid(us.3108, of_largest_polygon=TRUE), k=1), row.names=attr(us.poly_q, "region.id"), sym=TRUE)
us.poly_k <- union.nb(us.poly_q, us_k)
us.poly_k

You can check that this is a graph without subgraphs:

table(n.comp.nb(us.poly_k)$comp.id)

[Vickyxubba]

Thank you so much for the help! Really appreciate it!