前回はマルコフ連鎖モンテカルロ法にて熱中症搬送人員数をベイズ推定した.今回は空間相関を考慮した一般化線形モデル(ポアソンCARモデル)で熱中症搬送人員数をベイズ推定しようと試みた.結論から言うとうまく行っていない.途中でエラーが出てモデル構築に失敗する.その経過を記事とした.
spdepによる近接行列の生成
使用するRパッケージsf, NipponMap, spdepを呼び出す.
> library(NipponMap) > library(sf) Linking to GEOS 3.9.3, GDAL 3.5.2, PROJ 8.2.1; sf_use_s2() is TRUE > library(spdep) 要求されたパッケージ spData をロード中です To access larger datasets in this package, install the spDataLarge package with: `install.packages('spDataLarge', repos='https://nowosad.github.io/drat/', type='source')`
NipponMapのデータを読み込み,CRSを指定する.
> shp=system.file("shapes/jpn.shp", package = "NipponMap")[1] > pref=read_sf(shp) > st_crs(pref)=4326
EPSG4326では位置情報が緯度経度で与えられるため適切な距離計算が困難である.そのため投影変換により座標系を平面直角座標系第9系に変換する.
> pref=st_transform(pref, crs=6677)
ルーク型またはクイーン型の近接行列はspdepのpoly2nb関数で次のように生成する.しかしこの近接行列は後にnb2mat関数のところでエラーが出て使えなかった.
> nb=poly2nb(pref)
各都道府県の重心点の位置座標coordsを抽出し,最近隣4ゾーンで近接行列を定義する.knearneigh関数で探索して近接情報をknn形式のオブジェクトにまとめ,knn2nb関数でnb形式に変換する.
> coords=st_coordinates(st_centroid(pref)) Warning message: In st_centroid.sf(pref) : st_centroid assumes attributes are constant over geometries of x > knn=knearneigh(coords,4) > nb2=knn2nb(knn)
CARBasyesライブラリでベイズ推定
CARBayesライブラリを読み込む.
> library(CARBayes) 要求されたパッケージ MASS をロード中です 要求されたパッケージ Rcpp をロード中です Registered S3 method overwritten by 'GGally': method from +.gg ggplot2
ポアソンCARモデルの説明変数は日最高気温(x1),平均水蒸気圧(x2),ロジスティック関数を適用した月(x3),人口密度の対数(x4)とし,オフセット変数は65歳以上人口(p)の対数である.空間相関変数をLerouxモデルで与える場合,S.CARLeroux関数を用いて下記のように推定する.
> formula=y~offset(log(p))+x1+x2+x3+x4 > W=nb2mat(nb2, style = "B") > resLex=S.CARleroux(formula = formula, data=HeatStroke9, family = "poisson", W=W, burnin = 5000, n.sample = 20000)
しかし,結果は以下のようにエラーが出て失敗する.
Setting up the model. Error: W is not symmetric.
まとめ
ポアソンCARモデルにより空間相関を考慮した一般化線形モデルで熱中症搬送人員数をベイズ推定しようと試みたが,エラーのため完了できなかった.CARで用いる空間重み行列は対称行列でなければならないが,空間重み行列を対称化できない場合(例えば方向を考慮しなければならない場合など)CARの適用は不適切である(谷村晋『地理空間データ分析』p152).エラーの原因はCAR(条件付き自己回帰モデル)を適用しようとしたことにあるようだ.