SSブログ

hSDM.ZIP.iCARでうまくいかなかった例 [統計]

hSDMでIntrinsic CARをつかたZIPモデルをつかってみましたが、うまくいかなかったという例です。

データは、以前もつかったアラカシのデータです。コードは以下のとおり。

library(hSDM)
library(coda)

## Data
## http://cse.ffpri.affrc.go.jp/hiroki/stat/Qglauca.csv
qg <- read.csv("Qglauca.csv")

## Adjacent cells
adj <- c(2, 11,         # 1
         c(sapply(2:9, function(i) i + c(-1, 1, 10))),           # 2--9
         9, 20,         # 10

         c(sapply(seq(0, 170, 10),
                  function(j) j + c(1, 12, 21,
                                    c(sapply(2:9,
                                             function(i) i + c(0, 9, 11, 20))),
                                    10, 19, 30))), # 10--190
         181, 192,      # 191
         c(sapply(2:9, function(i) i + c(180, 189, 191))),     # 192--199
         190, 199)    # 200
num.adj <- c(2, rep(3, 8), 2,
             rep(c(3, rep(4, 8), 3), 18),
             2, rep(3, 8), 2)

qg$cell <- 1:nrow(qg)
pred <- data.frame(cell = qg$cell,
                   X = qg$X, Y = qg$Y)

fit.zip <- hSDM.ZIP.iCAR(counts = qg$N,
                         suitability = ~1,
                         abundance = ~1,
                         spatial.entity = qg$cell,
                         data = qg,
                         n.neighbors = num.adj,
                         neighbors = adj,
                         burnin = 10000,
                         mcmc = 10000,
                         thin = 10,
                         beta.start = 0,
                         gamma.start = 0,
                         Vrho.start = 10)
plot(fit.zip$mcmc)

しかし軌跡をみるとVrhoがどんどん発散していっているようです。うーむ。

Rplot.png


nice!(1)  コメント(0)  トラックバック(0) 
共通テーマ:日記・雑感

nice! 1

コメント 0

コメントを書く

お名前:
URL:
コメント:
画像認証:
下の画像に表示されている文字を入力してください。

Facebook コメント

トラックバック 0