如何从邮政编码创建一个邻接矩阵

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是最准确的。