如何解决Jags模型中的尺寸不匹配？

问题描述

我是贝叶斯分析领域的超级新手，我正在尝试使用经典捕获-捕获模型Mh2的示例进行练习

这是我的代码

nind <- dim(venados)[1] 
K <- 43 
ntraps <- 13
M <- 150 
nz <- M - nind 
Yaug <- array(0,dim = c(M,ntraps,K))
Yaug[1:nind,] <- venados
y <- apply(Yaug,c(1,3),sum) 
y[y > 1] <- 1

Bundle data
data1 <- list(y = y,nz = nz,nind = nind,K = K,sup = Buffer)

# Model JAGS 
sink("Mh2_jags.txt")
cat("
    model{

  # Priors
    p0 ~ dunif(0,1)
    mup <- log(p0/(1-p0))
    sigmap ~ dunif(0,10)
    taup <- 1/(sigmap*sigmap)
    psi ~ dunif(0,1)  

  # Likelihood
    for (i in 1:(nind+nz)) {   
    z[i] ~ dbern(psi)
    lp[i] ~ dnorm(mup,taup)
    logit(p[i]) <- lp[i]
    y[i] ~ dbin(mu[i],K)  
 } # i

    N <- sum(z[1:(nind+nz)])
    D <- N/sup*100
} # modelo
",fill = TRUE)
sink()

# Inicial values
inits <- function(){list(z = as.numeric(y >= 1),psi = 0.6,p0 = runif(1),sigmap = runif(1,0.7,1.2),lp = rnorm(M,-0.2))}

params1 <- c("p0","sigmap","psi","N","D")

# MCMC 
ni <- 10000; nt <- 1; nb <- 1000; nc <- 3

# JAGS and posteriors
fM2 <- jags(data1,inits,params1,"Mh2_jags.txt",n.chains = nc,n.thin = nt,n.iter = ni,n.burnin = nb)

我收到此错误消息

Processing function input....... 

Done. 
 
Compiling model graph
   Resolving undeclared variables
Deleting model

Error in jags.model(file = model.file,data = data,inits = inits,n.chains = n.chains,: 
  RUNTIME ERROR:
Compilation error on line 16.
Dimension mismatch in subset expression of y

我读到一些字母，如s和n必须更改。然而， 我不知道该怎么办。请提供意见。

非常感谢您

解决方法

问题是因为y是二维的，但是模型假定它是一维的。如果您假设次要调查是i.i.d. Bernoulli试用（每个会话都有K个试用）n，那么您只需要取y矩阵的行之和即可。假设是这种情况，那么您只需要在此脚本的顶部修改几行即可。

nind <- dim(venados)[1] 
K <- 43 
ntraps <- 13
M <- 150 
nz <- M - nind 
Yaug <- array(0,dim = c(M,ntraps,K))
Yaug[1:nind,] <- venados
y <- apply(Yaug,c(1,3),sum) 
y[y > 1] <- 1

# Take the rowSum
y_vector <- rowSums(y)

# Use y_vector instead of y
data1 <- list(y = y_vector,nz = nz,nind = nind,K = K,sup = Buffer)

相反，如果要包括观测过程的协变量（并且这些协变量因调查而异），则可以使用矩阵y并修改模型。

sink("Mh2_jags_Kloop.txt")
cat("
    model{

  # Priors
    p0 ~ dunif(0,1)
    mup <- log(p0/(1-p0))
    sigmap ~ dunif(0,10)
    taup <- 1/(sigmap*sigmap)
    psi ~ dunif(0,1)  

  # Likelihood
    for (i in 1:(nind+nz)) {   
    z[i] ~ dbern(psi)
    lp[i] ~ dnorm(mup,taup)
    logit(p[i]) <- lp[i]
      # Loop over K surveys
      for(j in 1:K){
        y[i,j] ~ dbern(p[i]*z[i])
    }  
 } # i

    N <- sum(z[1:(nind+nz)])
    D <- N/sup*100
} # modelo
",fill = TRUE)
sink()

最后，您没有指定模型中的mu。我认为您希望它是p，但是您还需要将潜在状态模型链接到观察状态模型（如果z=0则无法对该人进行采样。在这种情况下，您可以解释{{ 1}}作为psi个人在您网站上的概率。

nind+nz