为二项式数据编写分层JAGS模型时遇编译错误求助

给定二项式数据：

df <- data.frame(bin=c(1,1,1,1,1,0,0,0,0,0), group1 = c(1,1,0,0,0,0,0,0,0,0), group2 = c(0,0,1,1,1,0,0,0,0,0), fullgroup = c(1,1,2,2,2,0,0,0,0,0), precip = c(4,5,7,1,7,4,6,8,4,1), temp = c(1,5,8,2,4,7,3,9,2,4) )

尝试用JAGS构建分层logistic回归模型（公式为logit(p) = b0 + b1*precip + b2*temp），期望得到整体组（b0.overall、b1.overall、b2.overall）和分组（b0.group1、b1.group1等）的参数估计，但运行模型时出现编译错误：Compilation error on line 13. Attempt to redefine node b1[1]

错误原因 #

1. 笔误导致重复定义节点：模型中组参数的先验定义部分，若误将b2[j]的先验写成b1[j]，会重复定义b1[j]，触发该错误；同时原模型中b2[j]的先验错误使用了mu1.alpha/tau1.alpha，应改为对应b2的超参数mu2.alpha/tau2.alpha。
2. 未处理group=0的样本：JAGS数组索引从1开始，原模型未为group=0的样本指定参数，会导致未定义节点的问题，同时不符合"整体组"的建模需求。
3. 缺少似然定义：原模型仅计算了occurrence[i]，但未将响应变量y与occurrence[i]关联，模型无法完成推断逻辑。

解决方案 #

1. 修正模型笔误，确保每个参数的先验对应正确的超参数。
2. 为group=0的样本指定整体组参数（即超参数mu0.alpha、mu1.alpha、mu2.alpha），直接用这些值计算其线性预测值。
3. 补充响应变量y的似然定义，完善模型推断逻辑。
4. 修正初始值函数中rgamma的参数格式，添加组参数的初始值以优化收敛性。

修正后的JAGS模型 #

model {
  for(i in 1:length(y)) {
    # 区分group=0（整体组）和分组样本
    if(group[i] == 0) {
      psi[i] <- mu0.alpha + mu1.alpha*a[i] + mu2.alpha*b[i]
    } else {
      psi[i] <- b0[group[i]] + b1[group[i]]*a[i] + b2[group[i]]*b[i]
    }
    # 补充y的似然定义
    y[i] ~ dbern(1/(1+exp(-psi[i])))
  }
  
  for (j in 1:n.group) {
    b0[j] ~ dnorm(mu0.alpha, tau0.alpha)
    b1[j] ~ dnorm(mu1.alpha, tau1.alpha)
    b2[j] ~ dnorm(mu2.alpha, tau2.alpha)  # 修正为对应b2的超参数
  }
  
  # 超先验（对应整体组参数）
  mu0.alpha ~ dnorm(0, 0.00001)  # b0.overall
  tau0.alpha ~ dgamma(0.001, 0.001)
  mu1.alpha ~ dnorm(0, 0.00001)  # b1.overall
  tau1.alpha ~ dgamma(0.001, 0.001)
  mu2.alpha ~ dnorm(0, 0.00001)  # b2.overall
  tau2.alpha ~ dgamma(0.001, 0.001)
  
  # 可选：转换为标准差，方便结果解读
  sigma0.alpha <- 1/sqrt(tau0.alpha)
  sigma1.alpha <- 1/sqrt(tau1.alpha)
  sigma2.alpha <- 1/sqrt(tau2.alpha)
}

修正后的初始值函数 #

mod.inits <- function(){
  list(  
    mu0.alpha = rnorm(1, 0, 2),
    tau0.alpha = rgamma(n=1, shape=0.001, rate=0.001),
    mu1.alpha = rnorm(1, 0, 2),
    tau1.alpha = rgamma(n=1, shape=0.001, rate=0.001),
    mu2.alpha = rnorm(1, 0, 2),
    tau2.alpha = rgamma(n=1, shape=0.001, rate=0.001),
    b0 = rnorm(2, 0, 2),
    b1 = rnorm(2, 0, 2),
    b2 = rnorm(2, 0, 2)
  )
}

数据准备与模型运行 #

library(rjags)

data = list(
  y = df$bin,
  a = df$precip,
  b = df$temp,
  group = df$fullgroup,
  n.group = 2 
)

# 编译模型
ws <- jags.model(textConnection(model), data=data, n.chains = 3, n.adapt=30000)
# 额外燃烧期
update(ws, 10000)
# 抽取样本
samples <- coda.samples(ws, variable.names=c("mu0.alpha", "mu1.alpha", "mu2.alpha", "b0", "b1", "b2"), n.iter=50000)

内容的提问来源于stack exchange，提问作者Jennifer Jones