JAGS/R2Jags 中派生量的“收敛”

Question

更新：现在使用 Traceplot 示例

更新：现在有了新的跟踪图

我正在尝试调整 Outhwaite et. als 2018 代码以进行占用建模，并且有几个问题我似乎无法找到答案...

用于创建模型的代码

cat(
  "model{ 
  
  ### Model ###

        # State model
        for (i in 1:nsite){ 
          for (t in 1:nyear){   
            z[i,t] ~ dbern(psi[i,t]) 
            logit(psi[i,t])<- b[t] + u[i] 
          }}   
        
        # Observation model 
        for(j in 1:nvisit) {
          y[j] ~  dbern(Py[j]+0.0001)
          Py[j]<- z[Site[j],Year[j]]*p[j]      
          logit(p[j]) <- a[Year[j]] + c*logL[j]
        }
        
   ### Priors ###

      # State model priors
          for(t in 1:nyear){
            b[t] ~ dunif(-10,10)        # fixed year effect
          }                 
          
          for (i in 1:nsite) {
            u[i] ~ dnorm(0, tau.u)  # random site effect      
          } 
          
          tau.u <- 1/(sd.u * sd.u)  
          sd.u ~ dunif(0, 5)        # half-uniform hyperpriors
          
        
      # Observation model priors 
          for (t in 1:nyear) {
            a[t] ~ dnorm(mu.a, tau.a)   # random year effect           
          }
          
          mu.a ~ dnorm(0, 0.01)                         
          tau.a <- 1 / (sd.a * sd.a)    
          sd.a ~ dunif(0, 5)        # half-uniform hyperpriors             
          
          c ~ dunif(-10, 10)        # sampling effort effect
  
   
  ### Derived parameters ###
  # Finite sample occupancy - proportion of occupied sites
  for (t in 1:nyear) {  
    psi.fs[t] <- sum(z[1:nsite,t])/nsite
  } 
  
  #data# nyear, nsite, nvisit, y, logL, Site, Year
  
}", file="bmmodel.txt"
)

请注意，dbern(Py[j]+0.0001) 包含一个校正因子，因为 JAGS 不支持 dbern(0)。

我在一些工厂数据上运行模型，基本上只是尝试看看它是否像我预期的那样运行、收敛和表现。

问题 1（已回答）：我对 psi.fs[t] 的数量感兴趣。但由于模型在实际建模过程之后计算这个量，是否可以评估 psi.fs[t] 的收敛性？

使用R2JAGS运行模型的R代码

jagsrespsi<-jags(data.list, inits=test.inits,
     n.chains=2, n.iter=15000, n.thin=3,
     DIC=T, 
     model.file=paste0(modeltype,"model.txt"), parameters.to.save=c("psi.fs"))

问题 2：当我使用 traceplot(jagsrespsi) 绘制跟踪图时，似乎到处都是，但 jagsrespsi$BUGSoutput 的 Rhat 多年来一直是 1？ gelman.diag(as.mcmc(jagsrespsi)) 也表示收敛。监测 psi 也是如此！

我对这种模型行为感到非常惊讶，并怀疑有什么问题……但不知道该往哪里看

Answer 1

是的，您可以检查psi.ft[] 的收敛性，方法与检查模型参数的收敛性完全相同。这正是发生的情况，例如，在逻辑回归中，对于某些线性预测变量 z，拟合的响应概率计算为 exp(z)/(1 + exp(z))。

当你说跟踪图“到处都是”时，你是什么意思？这可能是好是坏。你能举个例子吗？一个“好”的跟踪图看起来像一条“肥大的毛毛虫”：从样本空间的所有区域中获取的连续样本，一个水平的毛球。尽管是为 SAS 编写的，this page 给出了一个合理的高级描述，描述了一个好的轨迹图是什么样的，以及不太理想的示例可能表明哪些问题。

响应您的编辑以包含跟踪图...

在我看来，这不是一个特别好的跟踪图：连续样本之间似乎存在一些负自相关。您计算出有效样本量 [ESS] 了吗？

但情节可能看起来有点奇怪，因为您的链条非常短，恕我直言。您可以使用 ESS 为估计概率的准确性提供非常粗略的近似值：二项式比例的最坏情况半宽 CI 为 +/-2 * sqrt(0.5*0.5/N)，其中 N 是样本大小（或本例中的 ESS） .因此，即使您的 MCMC 过程的效率为 1——因此 ESS 等于链长——那么您估计的准确度也仅为 +/-0.02。要将概率估计到小数点后 2 位（因此 CI 的半宽不超过 0.005），您需要 ESS 为 40,000。

在测试期间使用短链长度并没有错，但对于“生产”运行，我总是会使用远大于 2,500 的 chan 长度。 （而且我还会使用多个链，这样我就可以使用 Gelman-Rubin 统计来测试收敛性。）