JAGS - иерархическое сравнение моделей, не перепрыгивающее между моделями даже с псевдопримерами

Question

Я использую инфраструктуру иерархического моделирования, описанную Крушке, чтобы настроить сравнение между двумя моделями в JAGS. Идея в этой структуре состоит в том, чтобы запустить и сравнить несколько версий модели, указав каждую версию как один уровень категориальной переменной. Апостериорное распределение этой категориальной переменной может быть интерпретировано как относительная вероятность различных моделей.

В приведенном ниже коде я сравниваю две модели. Модели идентичны по форме. Каждый из них имеет один параметр, который необходимо оценить, mE, Как видно, модели отличаются по своим априорам. Оба приора распространяются в виде бета-версий, которые имеют режим 0,5. Тем не менее, предварительное распределение для модели 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

Хитрость заключается не в назначении фиксированной вероятности для модели, а в ее оценке (phi ниже) на основе единого априорного Затем вы хотите последующее распределение для phi поскольку это говорит вам о вероятности выбора модели 2 (т. е. "успех" означает m=1; Pr(модель 1) = 1-фи).

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

  m ~ dbern(phi)
  phi ~ dunif(0,1)

  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("phi")

Answer 2

Я не смог понять это, но я подумал, что любой, кто работает над этим, может по достоинству оценить следующий код, который будет воспроизводить проблему от начала до конца из R с помощью rjags (должен быть установлен 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)

Answer 3

Смотрите мой комментарий выше на ответ Марка С.

Этот ответ должен показать на примере, почему мы хотим сделать вывод по m, а не по phi.

Представьте, что у нас есть модель, заданная

data <- c(-1, 0, 1, .5, .1)

m~dbern(phi)
data[i] ~ m*dnorm(0, 1) + (1-m)*dnorm(100, 1)

Теперь очевидно, что истинное значение m равно 1. Но что мы знаем об истинном значении phi? Очевидно, более высокие значения phi более вероятны, но у нас нет достаточных доказательств, чтобы исключить более низкие значения phi. Например, phi=0.1 все еще имеет 10% -ную вероятность получения m=1; и phi = 0,5 все еще имеет 50% -ную вероятность получения m=1. Таким образом, у нас нет достаточных доказательств против довольно низких значений фи, хотя у нас есть железные доказательства того, что m=1. Мы хотим сделать вывод о m.