我在R2jags中看不到模型的结果--怎么了？

Question

我试图按照"Lahoz-Monfort JJ，Guillera-Arroita G，Tingley R (2015)统计方法“中提供的指令在R2jags中运行这个脚本。分子生态学资源，16,673-685。似乎一切正常，但我想不出要看结果的命令.有人能帮忙吗？

cat("model {
    # Priors
    psi ~ dunif(0,1)
    p11 ~ dunif(0,1)
    p10 ~ dunif(0,p10_max)
    
    # Likelihood 
    for (i in 1:S){
    z[i] ~ dbern(psi)
    p[i] <- z[i]*p11 + (1-z[i])*p10
    for (j in 1:K){
    Y[i,j] ~ dbern(p[i])
    }
    }
    } ",fill=TRUE)
sink()


Bayesian <- function(psi,p11,p10,S,K,nsims=100,doprint=TRUE,p10_max=0.05,
                                     ni=100000,nt=2,nc=1,nb=50000,myparallel=TRUE) {   
  psihat<-p11hat<-p10hat<-rep(nsims)
  modelSummaries<-list()
  for(ii in 1:nsims){
    if (doprint) cat("\r", ii, "of", nsims,"   ")    
    hh<-genSimData(psi,r11=0,p11,p10,S,K1=0,K2=K)
    # fit the model    
    jags.inits <-function()(list(psi=runif(1,0.05,0.95),p11=runif(1,p10_max,1),p10=runif(1,0,p10_max)))
    jags.data  <-list(Y=hh,S=S,K=K,p10_max=p10_max)
    jags.params<-c("psi","p11","p10")
    Thoropa_model<-jags(data = jags.data, inits = jags.inits, parameters.to.save= jags.params, 
                model.file= "Thoropa.txt", n.thin= nt, n.chains= nc, 
                n.iter= ni, n.burnin = nb, parallel=myparallel)  #, working.directory= getwd()
    # extract results (medians of the marginal posteriors)
    psihat[ii] <- model$summary["psi","50%"]
    p11hat[ii] <- model$summary["p11","50%"]
    p10hat[ii] <- model$summary["p10","50%"]    
    modelSummaries[[ii]]<-model$summary
  }
  if (doprint){
    printsummres(psihat,thename="estimated psi")
    printsummres(p11hat,thename="estimated p11")
    printsummres(p10hat,thename="estimated p10")
  }
  return(list(psihat=psihat,p11hat=p11hat,p10hat=p10hat,modelSummaries=modelSummaries))
  }

文件"Thoropa.txt“是一个存在/缺勤矩阵，如下所示：

PCR1    PCR2    PCR3    PCR4    PCR5    PCR6    PCR7    PCR8    PCR9    PCR10   PCR11   PCR12
1   0   0   0   0   0   0   0   0   0   0   0
0   0   0   1   0   0   0   1   0   0   0   0
0   0   0   0   0   0   1   0   0   0   0   0
1   1   1   1   1   1   1   1   1   1   1   1
0   0   0   0   0   0   0   0   0   0   0   0
1   0   1   0   1   1   1   1   1   1   1   1
0   0   1   0   0   0   1   0   0   0   0   0
1   0   1   0   0   0   0   0   1   0   0   0
0   0   0   0   0   0   0   0   0   0   0   0
1   0   0   0   1   1   0   0   0   1   0   0
1   1   0   1   0   1   0   1   0   0   1   0
1   1   0   0   0   0   0   1   0   0   0   0
0   0   0   0   0   0   0   0   0   0   0   0
1   1   1   1   0   1   1   1   1   0   1   1
1   1   1   1   1   1   1   1   1   1   1   1
0   0   1   0   1   0   0   0   0   0   1   1
1   1   1   1   1   1   1   1   1   1   1   1

以下是Limey的评论(谢谢！)我将脚本更改为：

sink("Thoropa2.txt")
cat("model {
    # Priors
    psi ~ dunif(0,1)
    p11 ~ dunif(0,1)
    p10 ~ dunif(0,p10_max)
    
    # Likelihood 
    for (i in 1:S){
    z[i] ~ dbern(psi)
    p[i] <- z[i]*p11 + (1-z[i])*p10
    for (j in 1:K){
    Y[i,j] ~ dbern(p[i])
    }
    }
    } ",fill=TRUE)
sink()

y=Thoropa# the detection/non detection table
S=nrow(y)
K=ncol(y)

psi ~ dunif(0, 1)
p11 ~ dunif(0, 1)
p10 ~ dunif(0, p10_max)
p10_max=0.05

jags.data<-list(y=y, S=S, K=K, p10_max=p10_max)

jags.inits <-function()(list(psi=runif(0,1),p11=runif(0,1),p10=runif(0,p10_max)))

jags.params<-c("psi","p11","p10")

Thoropa_model<-jags.parallel(data = jags.data, inits = jags.inits, parameters.to.save= jags.params, model.file= "Thoropa2.txt", n.chains= 4, n.thin= 10, n.iter = 100000, n.burnin=50000, jags.seed = 333) 

数据文件和以前一样。现在我收到了错误消息：

checkForRemoteErrors中的错误(Val)：4个节点产生错误；第一个错误:在边界外索引

有人能帮我识别我脚本中的错误吗？我不是专家，我在自学，很抱歉这是个愚蠢的问题.(也许数据格式有问题.？)

谢谢大家！

Answer 1

您的模型无法工作，因为您的R脚本中有一些语法错误。请注意，R语法与jags语法不同，即使您在de中运行jags，这些都是错误：

1. 符号"~“不用于R中的抽样。删除以下行：

(0，1) p11 ~ dunif(0，1) p10 ~ dunif(0，p10_max)

1. jags模型中的Y变量是大写的，所以必须纠正jags.data中的语法。

jags.data<-list(Y=y，S=S，K=K，K=K)

1. In jags.inits，a)函数的主体必须在大括号内，而不是括号内；b)函数runif有3个参数:n(要采样的值数)、min和max。正确的语法如下：

jags.inits <-function(){list(psi=runif(1,0,1)，p11=runif(1,0,1)，p10=runif(1,0，p10_max))}

修复这些错误，您的模型应该运行时没有错误。运行模型后，可以使用以下两个选项之一提取参数"psi“的中值：

   Thoropa_model$BUGSoutput$median$psi

   Thoropa_model$BUGSoutput$summary["psi","50%"]