2017-03-21 10 views
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 

しかし、そこには運。何か案は?

答えて

1

ここではループしないでijのインデックスを作成しています。この作業をするには、別のネストされたforループの中でそれを設定する必要があります。また、あなたの絶滅計算が間違っていました。

for(j in 2:T){ 
    for(i in 1:M){ 
    coln[i,j-1] <- ifelse(z[i,j]-z[i,j-1]==1,1,0) # colonized 
    ext[i,j-1] <- ifelse(z[i,j]-z[i,j-1]==-1,1,0) # went extinct 
} 
    tot.coln[j-1] <- sum(coln[,j-1]) # sum of colonized each survey 
    tot.ext[j-1] <- sum(ext[,j-1]) # sum of extinctions each survey 
    Nocc[j-1] <- sum(z[,j-1]) # total sites occupied each survey 
    coln.rate[j-1] <- tot.coln[j-1]/Nocc[j-1] 
    ext.rate[j-1] <- tot.ext[j-1]/Nocc[j-1] 
} 
+0

ご返信ありがとうございます。私はあなたのコードを「#導出量」行の後に置いたが、それはうまくいかなかった。私も 'for(i in 1:M)'ループを既存の '(j in 2:T)'ループ内に配置しようとしましたが、それも失敗しました。エラーが「ext.rateのサブセット式を評価できませんでした」と毎回。私はあなたの答えを間違って利用していますか? – tjr

+0

ああ。それは理にかなっている。 'tot.coln'、' tot.ext'、 'Nocc'、' coln.rate'、 'ext.rate'はすべて行列ではなくベクトルでなければなりません。コードを更新しました。このコードが動作しているかどうか教えてください。 –

+0

素晴らしい作品です、ありがとうございます! #Derived quantity linesの後にコードを配置しました。 – tjr