winbugs - JAGS - 即使使用伪先验，分层模型比较也不会在模型之间跳转

Question

我正在使用 Kruschke 描述的分层建模框架在 JAGS 中建立两个模型之间的比较。该框架的想法是通过将每个版本指定为分类变量的一个级别来运行和比较模型的多个版本。该分类变量的后验分布可以解释为各种模型的相对概率。

在下面的代码中，我比较了两个模型。这些模型在形式上是相同的。每个都有一个需要估计的参数，mE。可以看出，这些模型的先验有所不同。两个先验均以众数为 0.5 的 beta 分布形式进行分布。然而，模型 2 的先验分布更加集中。另请注意，我使用了伪先验，我希望它可以防止链卡在其中一个模型上。但无论如何，模型似乎都卡住了。

这是模型:

 model {

  m ~ dcat( mPriorProb[] )
  mPriorProb[1] <- .5
  mPriorProb[2] <- .5

  omegaM1[1] <- 0.5      #true prior
  omegaM1[2] <- 0.5      #psuedo prior 
  kappaM1[1] <- 3        #true prior for Model 1
  kappaM1[2] <- 5        #puedo prior for Model 1

  omegaM2[1] <- 0.5      #psuedo prior
  omegaM2[2] <- 0.5      #true prior
  kappaM2[1] <- 5        #puedo  prior for Model 2
  kappaM2[2] <- 10        #true prior for Model 2

  for ( s in 1:Nsubj ) {

    mE1[s] ~ dbeta(omegaM1[m]*(kappaM1[m]-2)+1 , (1-omegaM1[m])*(kappaM1[m]-2)+1 )
    mE2[s] ~ dbeta(omegaM2[m]*(kappaM2[m]-2)+1 , (1-omegaM2[m])*(kappaM2[m]-2)+1 )

    mE[s] <- equals(m,1)*mE1[s] + equals(m,2)*mE2[s]

    z[s] ~ dbin( mE[s] , N[s] )

  }
}

以下是相关数据的 R 代码:

dataList = list(
  z = c(81, 59, 36, 18, 28, 59, 63, 57, 42, 28, 47, 55, 38, 
        30, 22, 32, 31, 30, 32, 33, 32, 26, 13, 33, 30), 
  N = rep(96, 25),
  Nsubj = 25
)

当我运行此模型时，MCMC 的每次迭代都在 m = 1 处进行，并且从不跳转到 m = 2。我尝试了许多不同的先验和伪先验组合，但似乎找不到 MCMC 会考虑 m = 2 的组合。我什至尝试为模型 1 和 2 指定相同的先验和伪先验，但这没有帮助。在这种情况下，我预计 MCMC 会在模型之间相当频繁地跳转，大约花费一半的时间考虑一个模型，一半的时间考虑另一个模型。然而，JAGS 始终处于 m = 1 状态。我已经运行了长达 6000 次迭代的链，这对于像这样的简单模型来说应该足够长了。

如果有人对如何解决此问题有任何想法，我将不胜感激。

干杯，蒂姆

Answer 1

我还没能弄清楚这一点，但我认为从事这方面工作的其他人可能会欣赏下面的代码，它将使用 rjags 从 R 重现问题(必须安装 JAGS) .

请注意，由于此示例中只有两个竞争模型，因此我将 m ~ dcat() 更改为 m ~ dbern()，然后替换 m 与代码中其他地方的 m+1 。我希望这可以改善这种行为，但事实并非如此。另请注意，如果我们指定 m 的初始值，则无论我们选择哪个值，它都会停留在该值，因此 m 只是无法正确更新(而不是奇怪地被一个模型或另一个模型吸引)。这让我很头疼；为了马丁的眼睛，可能值得发布在 http://sourceforge.net/p/mcmc-jags/discussion/

library(rjags)
load.module('glm')

dataList = list(
  z = c(81, 59, 36, 18, 28, 59, 63, 57, 42, 28, 47, 55, 38, 
        30, 22, 32, 31, 30, 32, 33, 32, 26, 13, 33, 30), 
  N = rep(96, 25),
  Nsubj = 25
)

sink("mymodel.txt")
cat("model {

  m ~ dbern(.5)

  omegaM1[1] <- 0.5      #true prior
  omegaM1[2] <- 0.5      #psuedo prior 
  kappaM1[1] <- 3        #true prior for Model 1
  kappaM1[2] <- 5        #puedo prior for Model 1

  omegaM2[1] <- 0.5      #psuedo prior
  omegaM2[2] <- 0.5      #true prior
  kappaM2[1] <- 5        #puedo  prior for Model 2
  kappaM2[2] <- 10        #true prior for Model 2

    for ( s in 1:Nsubj ) {

    mE1[s] ~ dbeta(omegaM1[m+1]*(kappaM1[m+1]-2)+1 , (1-omegaM1[m+1])*(kappaM1[m+1]-2)+1 )
    mE2[s] ~ dbeta(omegaM2[m+1]*(kappaM2[m+1]-2)+1 , (1-omegaM2[m+1])*(kappaM2[m+1]-2)+1 )


    z[s] ~ dbin( (1-m)*mE1[s] + m*mE2[s] , N[s] )

    }
    }
    ", fill=TRUE)
sink()
inits <- function(){list(m=0)}

params <- c("m")

nc <- 1
n.adapt <-100
n.burn <- 200
n.iter <- 5000
thin <- 1
mymodel <- jags.model('mymodel.txt', data = dataList, inits=inits, n.chains=nc, n.adapt=n.adapt)
update(mymodel, n.burn)
mymodel_samples <- coda.samples(mymodel,params,n.iter=n.iter, thin=thin)
summary(mymodel_samples)