0

更新:现在使用 Traceplot 示例

迹线图

更新:现在有了新的跟踪图 正确的迹线图

我正在尝试适应Outhwaite 等。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) 包含一个校正因子,因为 dbern(0) 在 JAGS 中不受支持。

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

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

用于运行模型的 R 代码R2JAGS

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 也是如此!

我对这种模型行为感到非常惊讶,并且怀疑有什么问题......但不知道在哪里看

4

1 回答 1

1

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

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

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

对我来说,这看起来不是一个特别好的跟踪图:连续样本之间似乎存在一些负自相关。您是否计算过有效样本量 [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 统计来测试收敛性。)

于 2021-04-27T13:15:01.587 回答