2017-12-03 23 views
-1

私は、移植センターのための約208の郵便番号を持つ全米のための移植ケースに関するデータを持っています。私は緯度と経度の情報も持っています。私は、移植センターの空間的効果を見るために空間生存解析を試みています。 私はこのコードを持っており、近接性のためのadj.mtxの作成方法はわかりません。郵便番号で隣接行列を作成するにはどうすればいいですか

+1

私は質問に答えようとしましたが、あなた自身が求めているものがより明確であるように、質問自体に書式や文法の編集が必要です。 –

答えて

2

私は、実際の郵便番号が必ずしも物理的な隣接関係に対応していない数値隣接として有用であるとは思わない:ここでは

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

は私のデータがどのように見えるかです。あなたが求めているのは、郵便番号の近傍行列をどうやって取得するのかということだと思いますが、それを行うにはもっと情報が必要だと思います。たとえば、各郵便番号の中心(指定された緯度と経度)からの距離を計算してしきい値を設定できますが、一部の高密度エリアでは郵便番号の相対的なサイズによって距離が小さくなります。

geosphere package.このパッケージにはdistHaversineという関数呼び出しがあり、引数として2つの点とオプションの球面半径引数をとり、曲率を考慮して点間の距離を指定します地球とすべて。

これらのすべては、郵便番号の境界線でデータをダウンロードし、境界が交差する場所を調べる方がよいと言いました。私が見つけることができる最新のデータは、無料でhereです。それはかなり古いです。それにアクセスするには、ドロップダウンメニューをクリックして、郵便番号集計エリアを選択します。ファイルは約500MBです。

解凍されたデータには.shpファイルがあり、this postには形状データを扱うためのチュートリアルがあります。露骨にその男のポストからいくつかのものを取って、次のコードは、メモリにシェイプファイルを取得します。this答えと一緒に後

# 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) 

は、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) 

すべての郵便番号のためにそのコードを実行すると、非常に時間がかかり(およびおそらく場所での極端な違いに多少の誤差がスローされます)することができますが、アクセスにはあなたに興味がある唯一のものは、行は上のインデックスを作成しますことができます形状オブジェクト

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

自体は形状オブジェクト内のデータフレーム内にある実際の郵便番号を:例えば、最初の200のZIPコードの隣接関係を取得することによって行われます。

[email protected]$ZCTA5CE10 

をあなたがあなた自身のデータ内に存在する郵便番号のインデックスを取得するためにこれらを使用し、より迅速にその小さく、関連するサブセットに隣接行列を計算することができます:あなたがやってそれらにアクセスすることができます。このすべてでは、国勢調査データと独自のデータフレームを使用した簿記が必要になります。

購入可能な(または多分自由である)現在の郵便番号のジオメトリが多分存在する可能性がありますが、国勢調査データは利用可能で使いやすいものの良い組み合わせです。

2

Ryan Warnickに多分ポリゴンが必要だと私は同意します。しかし、maptools :: readShapePolyは使用しないでください。それは時代遅れの不完全な機能です。

使用あなたはライアンWarnickにより示唆されるように、あなた自身の隣接行列をロールバックしようとすることができますsfパッケージ(sf_read)から、またはrgdalからいずれかの機能(この場合はラスタパッケージを経由して容易になり)

library(raster) 
library(rgdal) 

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

しかし、そこには十分に確立された機能があります。特にspdepパッケージで。

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

あなたは、このようなspdep::nb2matとspdepなどの機能を持つnb :: sp2listw

P.S.を変換することができますgeosphereにいくつかの距離関数がありますが、distGeoが最も正確です。

+0

Upvoted。私はそれをするより良い方法があると考えていました。私はちょうどあなたがおそらくいくつかの外部データセットにアクセスする必要があることを渡りたいと思った。 –

関連する問題