Как я могу решить несоответствие размеров в модели Jags?

Я новичок в байесовском анализе и пытаюсь попрактиковаться на примере моделей Classic Capture-recapture: 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 должны быть изменены. Однако я не знаю, что делать. Пожалуйста, если вы могли бы дать совет.

Большое тебе спасибо


r bayesian jags
person SLV    schedule 13.08.2020    source источник

Ответы (1)


arrow_upward
0
arrow_downward

Проблема в том, что y является двумерным, но модель предполагает, что оно одномерное. Если вы предполагаете, что вторичные опросы являются i.i.d. испытаний Бернулли (а в каждом сеансе было 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, то этот человек не может быть отобран. В этом случае вы интерпретируете psi как вероятность того, что nind+nz человек находятся в вашем сайт.

# 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(p[i] * z[i],K)  
 } # i

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