0

我正在尝试按照“Lahoz-Monfort JJ, Guillera-Arroita G, Tingley R (2015) Statistical methods to account for false positive errors in environment DNA samples. Molecular Ecology Resources, 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 个节点产生错误；第一个错误：索引超出范围”

任何人都可以帮助确定我的脚本中的错误吗？我不是专家，我正在自学，如果这是一个愚蠢的问题，很抱歉......（也许数据格式有问题......？）

谢谢你们！

4

0 回答 0