如何從郵政編碼創建一個鄰接矩陣

Question

我有關於整個美國移植案例的數據，其中約208個郵政編碼用於移植中心。我也有他們的經緯度信息。我試圖進行空間生存分析，以查看移植中心的空間效應。 我有這個代碼，我不知道如何創建adj.mtx的接近度。

survregbayes(formula=Surv(time,cens)~age+sex + +frailtyprior("car",transplantcenter),data=d,survmodel="PO", + dist="loglogistic",mcmc=mcmc,prior=prior,Proximity=adj.mtx) 

這裏是我的數據是這樣的：

Answer 1

我不認爲實際的郵遞區號是作爲數字鄰接不一定對應於物理鄰接有用。我想你問的是如何獲得郵政編碼的鄰域矩陣，但要做到這一點，我相信你需要更多的信息。例如，您可以計算每個郵政編碼中心的距離（給定的緯度和經度）並將它們設爲閾值，但由於郵政編碼的相對大小，一些較高密度區域的距離會較小。

假設你想這樣做，你可以使用geosphere package.這個軟件包有一個函數調用distHaversine這需要作爲參數兩分，和一個可選的球面半徑參數，併爲您考慮到曲率的點之間的距離地球和一切。

所有這一切都說，更好的方法可能是下載郵政編碼邊界的數據並查看邊界相交的地方。我可以找到的最新數據（免費）是here;這是非常古老的。要訪問它，請單擊下拉菜單並選擇郵政編碼製表區。該文件大約500 MB。

解壓縮的數據有一個.shp文件，this post有一個很好的指導形狀數據的教程。公然採取了一些東西從那個傢伙的帖子，下面的代碼獲取形狀文件到內存：

# Install dependencies 
install.packages("rgeos") 
install.packages("maptools") 
library(rgeos) 
library(maptools) 
# Define the projection to be used 
crswgs84=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs") 
# Load the shapes 
postal.codes=readShapePoly("tl_2013_us_zcta510.shp",proj4string=crswgs84,verbose=TRUE) 

與this答案進行操作，請在rgeos包有一個功能叫做gTouches從多邊形構建鄰接矩陣邊界。

# Get adjacency matrix, returnDense=FALSE is to get a sparse matrix for memory purposes 
adj.mat <- gTouches(postal.codes, byid=TRUE, returnDense=FALSE) 

運行該代碼爲所有郵政編碼可以是非常耗時的（也可能引發一些錯誤，由於位置的極端差異），但只能訪問那些你有興趣，你可以做排索引上形狀對象。例如，得到的前200個郵政編碼的鄰接通過執行完成：

adj.mat.sub <- gTouches(postal.codes[1:200,], byid=TRUE, returnDense=TRUE) 

實際郵政編碼本身在形狀對象內的數據幀。您可以通過訪問它們：

postal.codes@data$ZCTA5CE10 

您可以使用這些來獲得指數出現在您自己的數據郵政編碼，以及更快速的計算上更小，相關子集的鄰接矩陣。所有這些都需要根據人口普查數據和您自己的數據框進行一些簿記。

可能有更多當前的郵政編碼幾何可供購買（或可能免費），但人口普查數據是可用且易於使用的良好組合。

Answer 2

我同意瑞恩Warnick，你可能需要多邊形。但是，您不應該使用maptools :: readShapePoly。這是一個過時和不完整的功能。

使用來自sf包（sf_read）或rgdal是函數（在這種情況下，通過光柵包更容易）

library(raster) 
library(rgdal) 

postal.codes <- shapefile("tl_2013_us_zcta510.shp") 

您可以嘗試推出自己的鄰接矩陣，由Ryan Warnick的建議，但是有很好的功能。值得注意的是在spdep包中。

library(spdep) 
nb <- poly2nb(postal.codes) 

您可以將nb與功能，如spdep::nb2mat和spdep :: sp2listw

附：在geosphere有幾個距離函數，但是distGeo是最準確的。