Depuración de JAGS y BUGS

Question

Estoy trabajando en un modelo bayesiano no muy elegante en R y JAGS. El objetivo es aislar los errores del codificador en una tarea de análisis de contenido. El código y la salida se dan a continuación.

La cuestión más importante es cómo depurar JAGS. (Asumo que el mismo consejo sería válido para BUGS también.) ¿Qué se supone que debo hacer con un error como "Valores iniciales inválidos" cuando hay casi una docena de valores iniciales diferentes?

Aquí está mi código R:

library(rjags)
library(R2jags)

#Load the data
toy_data <- read.csv("toy_data.csv")

#Rescale data
rescaled_data <- toy_data[,c(3:(2+K))]
for( k in c(1:K) ){
  col <- rescaled_data[,c(k)]
  rescaled_data[,c(k)] <- (col-min(col))/(max(col)-min(col))
}

codes <- as.matrix(rescaled_data)
doc_ids <- toy_data$docid
coder_ids <- toy_data$coderid

N <- nrow( toy_data ) #Number of document codings
D <- max(toy_data$docid) #Number of documents
I <- max(toy_data$coderid) #Number of coders
K <- dim(toy_data)[2]-2 #Number of attributes

#Package info for jags
jags.data <- list( "doc_ids", "coder_ids", "codes", "N", "D", "I", "K" )
jags.params <- c("z","mu","sigma","sigma_i","sigma_k","mu_dk","alpha_k","beta_k","alpha_i","beta_i")
jags.inits <- list(
  "z" <- matrix(rnorm(N*K),nrow=N,ncol=K),
  "mu" <- runif(1)*10,
  "sigma" <- rgamma(1,10),
  "sigma_i" <- rgamma(I,10),
  "sigma_k" <- rgamma(K,10),
  "mu_dk" <- as.matrix(rnorm(D*K),nrow=D,ncol=K),
  "alpha_k" <- rgamma(1,10),
  "alpha_i" <- rgamma(1,10),
  "alpha_i" <- rgamma(1,10),
  "beta_i" <- rgamma(1,10)
  )

#Fit the model
jagsfit <- jags(
  model.file="coder_model.txt",
  data=jags.data, 
  inits=jags.inits,
  jags.params,
  n.iter=5000,
)

Este es el modelo JAGS:

model {
  for( n in 1:N ){   #Loop over codings
    for( k in 1:K ){ #Loop over attributes
      #d <- doc_ids[n]    #Get document index
      #i <- coder_ids[n]  #Get code index
      codes[n,k] ~ dbern(p[n,k])
      logit(p[n,k]) <- z[n,k]
      z[n,k] ~ dnorm( mu_dk[doc_ids[n],k], sigma_k[k]*(1+sigma_i[coder_ids[n]]) )
    }
  }

  for( d in 1:D ){
    for( k in 1:K ){
      mu_dk[d,k] ~ dnorm( mu, sigma )
    }
  }

  for( k in 1:K ){   #Loop over attributes
    sigma_k[k] ~ dgamma( alpha_k, beta_k )
  }
  for( i in 1:I ){   #Loop over coders
    sigma_i[i] ~ dgamma( alpha_i, beta_i )
  }

  #Noninformative priors over alphas and betas
  mu ~ dnorm( 0, 10 )
  sigma ~ dgamma(10,8)
  alpha_k ~ dgamma(10,8)
  beta_k ~ dgamma(10,8)
  alpha_i ~ dgamma(10,8)
  beta_i ~ dgamma(10,8)
}

Estos son los datos:

"docid","coderid","Answer.1","Answer.2"
1,1,3,3
1,2,4,1
1,3,7,2
2,1,3,3
2,2,4,4
2,4,3,1
3,1,3,3
3,2,4,3
3,3,3,4
4,4,5,1
4,5,6,2
4,2,4,3
5,2,5,4
5,3,3,1
5,4,7,2
6,1,3,3
6,5,4,1
6,2,5,2

Y aquí está la salida de R:

Compiling model graph
   Resolving undeclared variables
   Allocating nodes
   Graph Size: 352

Error in jags.model(model.file, data = data, inits = inits, n.chains = n.chains,  : 
  Invalid initial values

Answer 1

He ejecutado su modelo con el paquete rjags. No he proporcionado ningún valor inicial ya que JAGS puede producirlos por ti. Usted puede ver el error a continuación

> m <- jags.model(file = "model.txt", n.chain = 1)
Compiling model graph
   Resolving undeclared variables
Deleting model

Error in jags.model(file = "model.txt", n.chain = 1) : RUNTIME ERROR:
Index out of range for node mu_dk

Ha especificado que la variable codes[n, k] en la probabilidad tiene una distribución Bernoulli como

codes[n, k] ~ dbern(p[n, k])

que debería tener valores de 0 o 1, pero en los datos tiene

> codes
      Answer.1  Answer.2
 [1,]     0.00 0.6666667
 [2,]     0.25 0.0000000
 [3,]     1.00 0.3333333
 [4,]     0.00 0.6666667
 [5,]     0.25 1.0000000
 [6,]     0.00 0.0000000
 [7,]     0.00 0.6666667
 [8,]     0.25 0.6666667
 [9,]     0.00 1.0000000
[10,]     0.50 0.0000000
[11,]     0.75 0.3333333
[12,]     0.25 0.6666667
[13,]     0.50 1.0000000
[14,]     0.00 0.0000000
[15,]     1.00 0.3333333
[16,]     0.00 0.6666667
[17,]     0.25 0.0000000
[18,]     0.50 0.3333333

¿Cómo pueden los códigos[n, k] tener una distribución Bernoulli?

Answer 2

Podría intentar ejecutar el mismo modelo en WinBugs u OpenBugs. Generalmente dan mensajes de error algo más detallados, por lo que también podrías obtener algo más útil en este caso concreto.