JAGS - hierarchical model comparison not jumping between models even with pseudopriors

Question

I'm using the hierarchical modelling framework described by Kruschke to set up a comparison between two models in JAGS. The idea in this framework is to run and compare multiple versions of a model, by specifying each version as one level of a categorical variable. The posterior distribution of this categorical variable then can be interpreted as the relative probability of the various models.

In the code below, I'm comparing two models. The models are identical in form. Each has a single parameter that needs to be estimated, mE. As can be seen, the models differ in their priors. Both priors are distributed as beta distributions that have a mode of 0.5. However, the prior distribution for model 2 is a much more concentrated. Note also that I've used pseudo priors that I had hoped would keep the chains from getting stuck on one of the models. But the model seems to get stuck anyway.

Here is the model:

 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] )

  }
}

Here is R code for the relevant data:

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
)

When I run this model, the MCMC spends every single iteration at m = 1, and never jumps over to m = 2. I've tried lots of different combinations of priors and pseudo priors, and can't seem to find a combination in which the MCMC will consider m = 2. I've even tried specifying identical priors and pseudo priors for models 1 and 2, but this was no help. In this situation, I would expect the MCMC to jump fairly frequently between models, spending about half the time considering one model, and half the time considering the other. However, JAGS still spent the whole time at m = 1. I've run chains as long as 6000 iterations, which should be more than long enough for a simple model like this.

I would really appreciate if anyone has any thoughts on how to resolve this issue.

Cheers, Tim

Answer 1

I haven't been able to figure this out, but I thought that anybody else who works on this might appreciate the following code, which will reproduce the problem start-to-finish from R with rjags (must have JAGS installed).

Note that since there are only two competing models in this example, I changed m ~ dcat() to m ~ dbern(), and then replaced m with m+1 everywhere else in the code. I hoped this might ameliorate the behavior, but it did not. Note also that if we specify the initial value for m, it stays stuck at that value regardless of which value we pick, so m just fails to get updated properly (instead of getting weirdly attracted to one model or the other). A head-scratcher for me; could be worth posting for Martyn's eyes at 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 2

The trick is not assigning a fixed probability for the model, but rather estimating it (phi below) based on a uniform prior. You then want the posterior distribution for phi as that tells you the probability of selecting model 2 (ie, a "success" means m=1; Pr(model 1) = 1-phi).

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 3

See my comment above on Mark S's answer.

This answer is to show by example why we want inference on m and not on phi.

Imagine we have a model given by

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

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

Now, it is obvious that the true value of m is 1. But what do we know about the true value of phi? Obviously higher values of phi are more likely, but we don't actually have good evidence to rule out lower values of phi. For example, phi=0.1 still has a 10% chance of yielding m=1; and phi=0.5 still has a 50% chance of yielding m=1. So we don't have good evidence against fairly low values of phi, even though we have ironclad evidence that m=1. We want inference on m.