WinBUGS、OpenBUGS中的类别变量

Question

我正在使用以下论文中的脚本(Zipkin，E.F.，Royle，J.A.，Dawson，D.K.，Bates，S.，2010 )。评估保护和管理行动的效果的多物种发生模型。生物保护143,479-484)来估计鸟类的占有率。我在检测估计中的一个变量(下面代码中的K循环)是Wind，它是一个分类变量，级别为1-6。我尝试在OpenBUGS中使用dcat函数，我希望它是一个统一的先验(beta(1,1))，但OpenBUGS失败并返回错误：

expected right parenthesis error pos 1344

当我删除这一行时，b3[i] ~ dcat(p[i,])#WIND Data不会发生。

任何关于如何正确指定dcat或如何为WinBUGS/OpenBUGS编写分类变量的建议都将不胜感激！

#model***************************************************************

occ.full <- function()
{

  #Define prior distributions for community-level model parameters
  omega ~ dunif(0,1)

  mean.u ~ dunif(0,1)
  mu.u <- log(mean.u) - log(1-mean.u)
  mean.v ~ dunif(0,1)
  mu.v <- log(mean.v) - log(1-mean.v)
  #mua1 ~ dnorm(0, 0.001)#WTR
  mua2 ~ dnorm(0, 0.001)#WBH
  mua3 ~ dnorm(0, 0.001)#GPH
  mua4 ~ dnorm(0, 0.001)#FHD
#   mua5 #OPEN this is percent data
  mua6 ~ dnorm(0, 0.001)#EdgeArea
  mua7 ~ dnorm(0, 0.001)#ForestArea
  mua8 ~ dnorm(0, 0.001)#Patch size
  mub1 ~ dnorm(0, 0.001)#OrdDate
  mub2 ~ dnorm(0, 0.001)#Start
  #mub3#WIND
  mub4 ~ dnorm(0, 0.001)#Temp

  tau.u ~ dgamma(0.1,0.1)
  tau.v ~ dgamma(0.1,0.1) 
  #tau.a1 ~ dgamma(0.1,0.1)
  tau.a2 ~ dgamma(0.1,0.1)
  tau.a3 ~ dgamma(0.1,0.1)
  tau.a4 ~ dgamma(0.1,0.1) 
#   tau.a5 #OPEN this is percent data
  tau.a6 ~ dgamma(0.1,0.1)#EDGE
  tau.a7 ~ dgamma(0.1,0.1)#Forest
  tau.a8 ~ dgamma(0.1,0.1) #Patch size
  tau.b1 ~ dgamma(0.1,0.1) 
  tau.b2 ~ dgamma(0.1,0.1)
  p ~ dbeta(1,1) #b3 is categorical
  tau.b4 ~ dgamma(0.1,0.1)

  for (i in 1:(n+nzeroes)) {
    #Create priors for species i from the community level prior distributions
    w[i] ~ dbern(omega)
    u[i] ~ dnorm(mu.u, tau.u)  
    v[i] ~ dnorm(mu.v, tau.v) 
    #a1[i] ~ dnorm(mua2, tau.a2)   
    a2[i] ~ dnorm(mua2, tau.a2)
    a3[i] ~ dnorm(mua3, tau.a3)
    a4[i] ~ dnorm(mua4, tau.a4)     
    a5[i] ~ dbeta(1, 1)#OPEN is percent data
    a6[i] ~ dnorm(mua6, tau.a6)#EdgeArea 
    a7[i] ~ dnorm(mua7, tau.a7)#ForestArea
    a8[i] ~ dnorm(mua8, tau.a8)#Patch size

    b1[i] ~ dnorm(mub1, tau.b1)    
    b2[i] ~ dnorm(mub2, tau.b2)
    b3[i] ~ dcat(p[i,])#WIND Data
    b4[i] ~ dnorm(mub4, tau.b4)

    #Create a loop to estimate the Z matrix (true occurrence for species i 
    #at point j.      
    for (j in 1:J) {
      logit(psi[j,i]) <- u[i] + a2[i]*WBH[j] + a3[i]*GPH[j] + a4[i]*FHD[j] + a5[i]*OPEN[j] + a6[i]*EdgeArea[j] + a7[i]*ForestArea[j] + a8[i]*Patch[j]
      mu.psi[j,i] <- psi[j,i]*w[i]
      Z[j,i] ~ dbern(mu.psi[j,i])

      #Create a loop to estimate detection for species i at point k during 
      #sampling period k.      
      for (k in 1:K[j]) {  
        logit(p[j,k,i]) <- v[i] + b1[i]*OrdDate[j,k] + b2[i]*Start[j,k] + b4[i]*Temp[j,k] + b3[i]*Wind[j,k] 
        mu.p[j,k,i] <- p[j,k,i]*Z[j,i]
        X[j,k,i] ~ dbern(mu.p[j,k,i])
      }#K   
    }#J
  }#N

  #Sum all species observed (n) and unobserved species (n0) to find the 
  #total estimated richness
  n0 <- sum(w[(n+1):(n+nzeroes)])
  N <- n + n0

  #Create a loop to determine point level richness estimates for the 
  #whole community and for subsets or assemblages of interest.
  for(j in 1:J){
    Nsite[j]<- inprod(Z[j,1:(n+nzeroes)],w[1:(n+nzeroes)])
    Nfor[j]<- inprod(Z[j,1:n],h.for[1:n])
    Nnon[j]<- inprod(Z[j,1:n],h.non[1:n])
    Nint[j]<- inprod(Z[j,1:n],h2.int[1:n])
    Nedge[j]<- inprod(Z[j,1:n],h2.edge[1:n])
    #Nneo[j]<- inprod(Z[j,1:n],m.neo[1:n])
  }
}

Answer 1

我不确定dcat发行版本身，但您的错误可能是因为您有：

b3[i] ~ dcat(p[i,])#WIND Data

但是，当您在模型中指定p时，它是一个三维数组p[j,k,i].，我认为您需要：

b3[i] ~ dcat(p[,,i])#WIND Data

请注意，i是p的最后一个索引，两个逗号。希望他能帮上忙。