r - JAGS を使用して、2 項および 20 サンプル/プロットで林冠被覆の中央値を推定する

Question

13 のグループ (植生の有無にかかわらず 20 からのポイント数として記録) と合計 223 のサンプルについて、中央値と、5、25、75、および 95 パーセンタイルのキャノピー カバーを推定したいと思います。以前、ベータ版配布を想定してこれを投稿しましたが、それは正しくありませんでした。これは提出期限が過ぎた原稿のためのもので、これが最後に欠落している部分です。誰かが私を完成させるのを手伝ってくれたら本当にありがたいです（コードが機能するまで）。私はそれに近いと思いますが、微調整が必​​要です-私は思います。

（私は編集して2つの反対票を修正しましたが、何が明確でないのかわかりません）。

大変感謝します！

以下は、私のモデル ステートメント、R コード、およびデータです。そのまま、私が得るエラーは

model(model.file, data = data, inits = init.values, n.chains = n.chains,  : 
  RUNTIME ERROR:
Compilation error on line 16.
Subset out of range: re[14]

ただし、以下のスペースを削除したことに注意してください。エラーは可能性のステートメントを参照しています。

model{
# priors
for (i in 1:13){
 alpha[i] ~ dunif(0, 1)
 re[i] ~ dnorm(0, 0.001)
}
#likelihood
for (i in 1:223) {
   canopy[i] ~ dbin(p[i], 20)
   logit(p[i]) <- alpha[site[i]] + re[i] 
}   
 median <- 1/(1+exp(alpha[site[i]]))
t4est1_100  <- step(median[1]-median[4])
t5est1_10 <- step(median[3]-median[4])
t6est10_100 <- step(median[2]-median[3])

}

R コード:

cover <- read.csv("f:\\brazil\\canopy2.csv", header=T)
library(R2jags)
library(rjags)
setwd("f://brazil")
site <- frag$site
canopy <- frag$canopy*20
N <- length(frag$site)

jags.data <- list("site", "canopy")
jags.params <- c("median", "test100MF","test100MT","test100fc","test100fa", 
"test100gv","test100hm","test100mc", "test100ca","test100ct", "test10MF",
"test10MT", "test10fc","test10fa", "test10gv", "test10hm", "test10mc", "test10ca", 
"test10ct", "test1MF", "test1MT", "test1fc",  "test1fa",  "test1gv", "test1hm", 
"test1mc", "test1ca", "test1ct", "t1est1_con","t2est10_con","t3est100_con",
"t4est1_100","t5est1_10","t6est10_100")
#inits1 <- list(a=0, sd=0)
#inits2 <- list(a=100, sd=50)
#jags.inits <- list(inits1, inits2)

jags.inits <- function() {
  list(alpha = 0, re=0)}

jagsfit2 <- jags(data=jags.data, inits=jags.inits, jags.params,
n.iter=1000000, n.burnin=20000, model.file="fragmodelbinom.txt")
my.coda <- as.mcmc(jagsfit2)
summary(my.coda, quantiles=c(0.05, 0.25,0.5,0.75, 0.95))
print(jagsfit2, digits=3)

データ：

structure(list(site = c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 
4L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 
5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 6L, 6L, 6L, 6L, 
6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 
7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 
8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 9L, 9L, 9L, 9L, 10L, 10L, 
10L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 12L, 12L, 12L, 12L, 
13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L, 13L
), canopy = c(0, 0.05, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 
0.6, 0.6, 0.65, 0.65, 0.7, 0.7, 0.7, 0.7, 0.7, 0.7, 0.75, 0.75, 
0.8, 0.8, 0.85, 0.85, 0.85, 0.85, 0.85, 0.85, 0.9, 0.9, 0.9, 
0.9, 0.95, 0.95, 0.95, 0.95, 0.95, 1, 1, 1, 1, 1, 1, 0.05, 0.2, 
0.25, 0.4, 0.4, 0.5, 0.6, 0.6, 0.65, 0.65, 0.75, 0.75, 0.75, 
0.8, 0.8, 0.8, 0.8, 0.85, 0.85, 0.85, 0.9, 0.9, 0.95, 0.95, 0.95, 
0.95, 1, 1, 1, 1, 1, 0, 0.2, 0.25, 0.3, 0.35, 0.4, 0.4, 0.45, 
0.45, 0.5, 0.5, 0.55, 0.6, 0.7, 0.7, 0.7, 0.8, 0.85, 0.9, 0.9, 
0.9, 0.95, 0.95, 0.95, 0.95, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 
0.1, 0.4, 0.4, 0.45, 0.5, 0.55, 0.55, 0.7, 0.7, 0.75, 0.8, 0.8, 
0.8, 0.9, 1, 1, 0.15, 0.2, 0.25, 0.25, 0.35, 0.5, 0.5, 0.55, 
0.65, 0.7, 0.7, 0.75, 0.8, 0.85, 0.85, 0.9, 0.9, 0.95, 0.95, 
1, 1, 1, 1, 0.05, 0.4, 0.6, 0.65, 0.65, 0.65, 0.7, 0.85, 0.95, 
1, 1, 1, 0.35, 0.4, 0.4, 0.5, 0.5, 0.55, 0.65, 0.65, 0.75, 0.75, 
0.8, 0.85, 0.9, 0.9, 1, 1, 1, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7, 
0.8, 0.8, 0.8, 0.85, 0.95, 0.95, 1, 1, 1, 1, 0.8, 0.85, 1, 1, 
1, 1, 1, 0, 0.05, 0.1, 0.15, 0.5, 0.6, 0.6, 0.75, 0.1, 0.35, 
0.6, 1, 0.4, 0.5, 0.55, 0.65, 0.65, 0.8, 0.9, 0.9, 0.9, 0.9, 
0.95, 0.95, 1)), .Names = c("site", "canopy"), class = "data.frame", 
row.names = c(NA, -227L)) 

Answer 1

JAGS では、dbin成功確率と試行回数を取り、成功回数を示す整数を返します。モデルでは、結果がすべて 0 と 1 の間にある値であることを指定していますがdbin、canopyこれは正しくありません。

なぜ 20 回の試行を指定するのですか? 種の有無を評価した 20 点はありますか? その場合、おそらくデータ モデルを次のように指定する必要があります。

y[i] ~ dbin(cover[i], 20)

ここで、y[i]はその場所に存在する種の数を示す整数i(これらはあなたのデータです) であり、coverは推定しようとしているものです。つまり、観察として記録された種の割合のカバーですi)。( Fukaya et al. (2010)は同様の問題のコードを提供しており、調査区画でのフジツボの発生確率を推定しています。一読の価値があります。)

これは、あなたの問題 (または少なくとも私の理解) を表すシミュレートされたデータを使用した例です。これはコンパイルされ、JAGS がかなり正確なカバー値を推定します。

M <- function() {
  for (j in 1:13){
    alpha[j] ~ dnorm(0, 0.0001)
  }

  for (i in 1:length(y)) {
    y[i] ~ dbin(p[i], 20)
    logit(p[i]) <- alpha[site[i]] + eps[i]
    eps[i] ~ dnorm(0, sd^-2)
  }   
  sd ~ dunif(0, 100)
}

はロジット スケール上にあり、0 と 1 に制限されていないためalpha、uniform(0, 1) の事前ではなく、漠然とした法線の事前を に使用したことに注意してください。わかりやすくするために に名前を変更し、移動しました観測レベルのループに。alphareeps

site <- c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
          1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
          1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 
          2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 
          3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 
          3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 
          4, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 
          5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 
          7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 8, 8, 8, 8, 8, 8, 8, 8, 8, 
          8, 8, 8, 8, 8, 8, 8, 9, 9, 9, 9, 10, 10, 10, 11, 11, 11, 11, 
          11, 11, 11, 11, 12, 12, 12, 12, 13, 13, 13, 13, 13, 13, 13, 13, 
          13, 13, 13, 13, 13)

set.seed(1)
alpha <- rnorm(13) # 13 independent random site intercepts
sd <- 0.123        # a constant variance for obs around alphas
logit.p <- rnorm(length(site), alpha[site], sd) # logit probs
y <- sapply(plogis(logit.p), rbinom, n=1, size=20)

シミュレートしたデータを示す 2 つのプロット:

plot(plogis(logit.p) ~ site, pch=20, las=1, ylab='p', ylim=c(0, 1))

plot(y ~ jitter(site), pch=20, las=1, xlab='site')

fit <- jags(list(y=y, site=site), NULL, c('sd', 'alpha', 'p'), M, 3, 1000)

pJAGS の推定値との比較:

plot(plogis(logit.p), fit$BUGSoutput$summary[
  grep('^p\\[', row.names(fit$BUGSoutput$summary)), '50%'], pch=20,
  xlab='True p', ylab='Estimated p (median and 95% CrI)', las=1, ylim=c(0, 1),
  panel.first=abline(0, 1, lwd=3, col='light gray'))
segments(plogis(logit.p), 
         fit$BUGSoutput$summary[
           grep('^p\\[', row.names(fit$BUGSoutput$summary)), '2.5%'], 
         y1=fit$BUGSoutput$summary[
           grep('^p\\[', row.names(fit$BUGSoutput$summary)), '97.5%'])

そしてアルファ：

plot(alpha, fit$BUGSoutput$summary[
  grep('^alpha\\[', row.names(fit$BUGSoutput$summary)), '50%'], pch=20,
  xlab='True alpha', ylab='Estimated alpha (median +95% CrI)', las=1, 
  ylim=c(-1.5, 2.5), panel.first=abline(0, 1, lwd=3, col='light gray'))
segments(alpha, 
         fit$BUGSoutput$summary[
           grep('^alpha\\[', row.names(fit$BUGSoutput$summary)), '2.5%'], 
         y1=fit$BUGSoutput$summary[
           grep('^alpha\\[', row.names(fit$BUGSoutput$summary)), '97.5%'])

要約 (分位数など) は保存されたサンプルから簡単に抽出でき、これらの一部はすでに JAGS 要約 (つまりfit$BUGSoutput$summary) に表示されています。あなたの場合、推定の分位点alpha(つまり、サイト カバーの平均) におそらく関心があるでしょう。これらにアクセスする別の方法は、 によってまだ計算されていない分位数が必要な場合に役立ちますR2jags。次の方法があります。

apply(fit$BUGSoutput$sims.list$alpha, 2, quantile, c(0.2, 0.5, 0.8))

#           [,1]       [,2]       [,3]     [,4]      [,5]       [,6]      [,7]
# 20% -0.6753085 0.08214699 -0.8994924 1.213293 0.2411191 -0.9441959 0.4038927
# 50% -0.6151128 0.15716925 -0.8299434 1.310542 0.3214097 -0.8205142 0.4909790
# 80% -0.5461859 0.23436643 -0.7632734 1.405498 0.3994087 -0.6982591 0.5887423
#          [,8]      [,9]       [,10]    [,11]     [,12]      [,13]
# 20% 0.6052447 0.6722630 -0.46050504 1.513086 0.2352861 -0.6516170
# 50% 0.6971180 0.8914973 -0.24458871 1.677613 0.4217723 -0.5443512
# 80% 0.7858756 1.1142837 -0.01018129 1.839783 0.6183198 -0.4309478

alphaただし、共通の分布に属するアルファの母集団から抽出される、階層的なアプローチを検討することをお勧めします。これには次のようなものが含まれます。

M <- function() {
  for (j in 1:13){
    alpha[j] ~ dnorm(mean.alpha, sd.alpha^-2)
  }

  for (i in 1:length(y)) {
    y[i] ~ dbin(p[i], 20)
    logit(p[i]) <- alpha[site[i]] + eps[i]
    eps[i] ~ dnorm(0, sd^-2)
  }   
  sd ~ dunif(0, 100)
  mean.alpha ~ dnorm(0, 0.0001)
  sd.alpha ~ dunif(0, 100)
}

これにより、グループ平均でのサイト間の分散について何かがわかります (つまり、 を検査することによってsd.alpha)。また、データが少なすぎる可能性がありますが、サイト間で異なる (つまり、サイト内の観測間の分散) を許可することを検討することもできsdます (とにかく、単一の定数を持つことは完全に賢明sdです... システムによって異なります)。 .

Answer 2

コードを読むと、明らかな間違いが 2 つあります。モデル定義には次のものがあります。

    for (i in 1:223) {
     canopy[i] ~ dbin(p[i], 20)
     logit(p[i]) <- alpha[site[i]] + re[i] 
    } 

re定義したように、パラメーターには 13 個のコンポーネントしかないため、上記は正しくありません。re[site[i]]代わりにすべきだと思います。これが、発生しているエラーの原因です。

2 番目は行median <- 1/(1+exp(alpha[site[i]]))です。これは for ループの外側にありますが、に依存しているため、内側にある必要があると思いますi。