OpenBUGS:无法编译分层模型

Question

我在OpenBUGs中使用以下数据实现了一个模型：

​

​

数据由10项独立研究组成-2种不同药物(Med =0或1)的5项试验，然后是该试验(感染)的全部感染参与者-以及试验中的总参与者。

bugsData(list(y = df$Infected, n = 10, J = 5 , t = df$Total, x = df$Med) , file="Bugs_Data.txt")

我的原始模型可以写成：

Y[i]|mu[i],t[i] ~ Bin(mu[i], t[i]), i = 1,.....10
logit(mu[i]) = Beta0 + Beta1*x[i]

我可以像这样在OpenBUGS中实现它：

model {
## LIKELIHOOD
  for (i in 1:n) {
    y[i] ~ dbin(mu[i], t[i])
    logit(mu[i]) <- Beta0 + Beta1*(x[i])
  }
  ## NORMAL PRIORS 
  Beta0 ~ dnorm(0, 0.0001)
  Beta1 ~ dnorm(0, 0.0001)
}

然而，我想实现一个分层模型，其中s(i)表示来自列试验的试验，其中进行了观察I：

Y[i]|mu[i],t[i] ~ Bin(mu[i], t[i]), i = 1,.....10
logit(mu[i]) = Beta0,s(i) + Beta1*x[i]
Beta0,j ~ Normal(mu_bo, sigma2_b0), j = 1,....5

我已经尽我所能地尝试了一个模型，但到目前为止还没有成功。

model {
## LIKELIHOOD
  for (j in 1:J) {
    for (i in 1:n) {
      y[i] ~ dbin(mu[i], t[i])
      logit(mu[i]) <- Beta0[j] + Beta1*(x[i])
    Beta0[j] ~ dnorm(mu_b0, tau_b0)
    }
  }
  
  ## PRIORS 
  mu_b0 ~ dnorm(0, 0.0001)
  tau_b0 ~ dgamma(0.0001, 0.0001)
  sigma2_b0 <- 1 / tau_b0
  Beta1 ~ dnorm(0, 0.0001)
}

当前模型不能按原样使用我的数据进行编译。

Answer 1

您需要提供一个长度为n的向量，用于索引i所属的组数据点。

model {
## LIKELIHOOD
for (i in 1:n) {
  y[i] ~ dbin(mu[i], t[i])
  logit(mu[i]) <- Beta0[group_vec[i]] + Beta1*(x[i])
}
## PRIORS 
mu_b0 ~ dnorm(0, 0.0001)
tau_b0 ~ dgamma(0.0001, 0.0001)
sigma2_b0 <- 1 / tau_b0
Beta1 ~ dnorm(0, 0.0001)
for(j in 1:J){
  Beta0[j] ~ dnorm(mu_b0, tau_b0)
}
}

group_vec的长度为n，它采用i所属的组数据点的整数。例如，如果我们有21个数据点和3个组。

groups <- factor(rep(c("A", "B", "C"), each = 7))

然后，我们可以将group_vec构造为

group_vec <- as.numeric(groups)
[1] 1 1 1 1 1 1 1 2 2 2 2 2 2 2 3 3 3 3 3 3 3

这被称为嵌套索引，这是指定分层模型的一种非常常见的方法，您只需将其与数据一起提供即可。在你的例子中，你已经有了一个整数形式的试验。因此，您可以将group_vec作为数据中的试验列，也可以将模型中的group_vec重命名为Trial。