Блог

У меня есть фрейм данных, такой как ;

 Orders  Time  NbBones
1      12    21
1      34    43
2      32    0
3      21    0
2      10    0
 

И модель джагса пуассона:

 Poisson_model<-"
model 
{
  for(i in 1 : N) 
  {
    NbBones[i] ~ dpois(lambda[i])
    lambda[i] <-  alpha   Coeforder[Orders[i]]    Coeftime * Time[i] 
  }
  Coeforder[1] <- 0
  Coeforder[2] ~ dnorm(0,1)
  Coeforder[3] ~ dnorm(0,1) 
  Coeftime ~ dnorm(0,1) 
  alpha ~ dexp(1)  
}
"

data4jags <- list(Time = data$Time, NbBones = data$NbBones, Orders=data$Orders, N=length(data$Orders))


m1 <- jags.model(textConnection(Poisson_model), data = data4jags, n.chains = 3)
 

Но когда я запускаю код в R, я получаю следующее сообщение об ошибке :

 Error in jags.model(textConnection(Model1), data = data4jags, n.chains = 3) : 
  Error in node NbBones[4]
Node inconsistent with parents
 

Кто-нибудь понимает, в чем проблема в коде ? Большое спасибо за ваше время

Для Бена Болкера

Ты имеешь в виду такой код ?

 model1<-"
model 
{
  for(i in 1 : N) 
  {
    NbEVES_branche[i] ~ dpois(lambda[i])
    partial[i] <- alpha   Coeftime * Time[i]
    if (Orders[i]>1) {
    lambda[i] <-  partial[i]   Coeforder[Orders[i-1]]
}
Coeforder[1] ~ dnorm(0,1)
Coeforder[2] ~ dnorm(0,1) 
Coeftime ~ dnorm(0,1) 
alpha ~ dexp(1)  

  }
  "
 

Две вещи бросаются мне в глаза. (1) ваша модель допускает отрицательные значения для среднего значения Пуассона (хотя я ожидаю, что это даст вам другую ошибку); (2) Я был бы удивлен, если бы вам разрешили использовать массив ( Coeforder ), где первый элемент является логическим узлом (0), в то время как остальные являются стохастическими узлами. В конце концов это сработало для меня:

 dd <- read.table(header=TRUE, text = "
Orders  Time  NbBones
1      12    21
1      34    43
2      32    0
3      21    0
2      10    0
")


ddj <- c(as.list(dd), N=nrow(dd))

model1<-"model {
  for(i in 1 : N) {
    NbBones[i] ~ dpois(lambda[i])
    partial[i] <- alpha   Coeftime * Time[i]
    lambda[i] <- ifelse(Orders[i] > 1,
         partial[i]   Coeforder[max(1,Orders[i]-1)],
         partial[i])
    }
    Coeforder[1] ~ dnorm(0,1)
    Coeforder[2] ~ dnorm(0,1) 
    Coeftime ~ dnorm(0,1)
    alpha ~ dexp(1)
}"

writeLines(model1, "tmp.jags")
library(rjags)
m1 <- jags.model("tmp.jags", ddj, n.chains = 3)
 

Обратите внимание, что у JAGS нет if/else заявления. Это max(1,...) предложение необходимо, потому ifelse() что оператор всегда оценивает оба своих потенциальных результата, даже если они не будут использоваться (поэтому мы должны быть осторожны, чтобы не ссылаться на Coeforder[0] них ).

Более компактно, вы можете настроить матрицу модели X вне JAGS и использовать структуру, например lambda[i] <- exp(inprod(X[i,], beta)) , где beta находится ваш полный вектор коэффициентов (это имитирует структуру стандартного пуассоновского GLM, включая компонент логарифмической связи/экспоненциальной обратной связи).