1
私は比較的JAGSが新しく、RパッケージのjagsUIを通して実行しています。私は占有モデルを構築していますが、結果を要約したいと思います。だから私は0と1の行列いますJAGSポスト計算とifelse/step
mat1 <- matrix(rbinom(10*10,1,.5),10,10)
y=mat1
私は以下のモデルを介して実行したい:このモデルは、「#派生量」の後に作品を除いて
# Bundle data and summarize data bundle
str(win.data <- list(y = mat1, M = nrow(mat1), T = ncol(mat1)))
# Specify model in BUGS language
sink("model.txt")
cat("
model {
# Priors
psi0 ~ dunif(0, 1)
p ~ dunif(0, 1)
for(t in 1:(T-1)){
rho[t] ~ dunif(-1,1)
}
beta0 ~ dnorm(0, 0.1)
# Likelihood
for (i in 1:M) { # Loop over sites
z[i,1] ~ dbern(psi0) # State model
y[i,1] ~ dbern(z[i,1]*p)
for (j in 2:T) { # Loop over replicate surveys
logit(psi[i,j])<- beta0 + rho[j-1]*z[i,j-1]
z[i,j] ~ dbern(psi[i,j])
y[i,j] ~ dbern(z[i,j]*p) # Observation model
}
}
# Derived quantities
coln[i,j] <- ifelse(z[i,j]-z[i,j-1]==1,1,0) # colonized
ext[i,j] <- ifelse(z[i,j-1]-z[i,j]==1,1,0) # went extinct
tot.coln[,j] <- sum(coln[,j]) # sum of colonized each survey
tot.ext[,j] <- sum(ext[,j]) # sum of extinctions each survey
Nocc[,j] <- sum(z[,j]) # total sites occupied each survey
coln.rate[,j] <- tot.coln[,j]/Nocc[,j]
ext.rate[,j] <- tot.ext[,j]/Nocc[,j]
}
",fill = TRUE)
sink()
# Initial values
zst <- apply(y, 1, max, na.rm=TRUE) # Avoid data/model/inits conflict
y<- as.matrix(y)
zst<- y
inits <- function(){list(z = zst)}
# Parameters monitored
params <- c("psi0", "p", "beta0", "coln.rate", "ext.rate")
# MCMC settings
ni <- 2000 ; nt <- 1 ; nb <- 1000 ; nc <- 3
# Call JAGS and summarize posteriors
library(jagsUI)
fm <- jags(win.data, inits, params, "model.txt", n.chains = nc,
n.thin = nt, n.iter = ni, n.burnin = nb)
print(fm, dig = 3)
を実行します。基本的には、各調査で0から1まで、1から0までの変化率を計算したいと思います。それがうまくいかない理由についての私の考えのカップル。 1)z [i、j]は実際には0と1ではありません。 2)計算は派生量にはならない。 3)JAGSマニュアルのifelseが私の考えをしていない。
私もで最初の2行を置換した後派生量を「ステップ」機能を使用してみました:
coln[i,j] <- step(z[i,j]-z[i,j-1]-0.5) # colonized
ext[i,j] <- step(z[i,j-1]-z[i,j]-0.5) # went extinct
しかし、そこには運。何か案は?
ご返信ありがとうございます。私はあなたのコードを「#導出量」行の後に置いたが、それはうまくいかなかった。私も 'for(i in 1:M)'ループを既存の '(j in 2:T)'ループ内に配置しようとしましたが、それも失敗しました。エラーが「ext.rateのサブセット式を評価できませんでした」と毎回。私はあなたの答えを間違って利用していますか? – tjr
ああ。それは理にかなっている。 'tot.coln'、' tot.ext'、 'Nocc'、' coln.rate'、 'ext.rate'はすべて行列ではなくベクトルでなければなりません。コードを更新しました。このコードが動作しているかどうか教えてください。 –
素晴らしい作品です、ありがとうございます! #Derived quantity linesの後にコードを配置しました。 – tjr