transdimensional model switch index is not jumping at all

Help
2013-10-09
2013-10-30
  • Hi folks,
    I'm trying to implement a model switch approach as described by Lodewyckx (2011) and Kruschke in Chapter 8 of his textbook.
    The basic idea is to set up two models in one supermodel and sample a model switching index to get an estimation for the bayes factor. My problem is, that the model index is not jumping at all in 100000 iterations (JAGS 3.4.0). Neither the bisection method of Lodewyckx nor setting up pseudo priors gave any help.

    Do you have any suggestions how I could get the index to switching? Or another possibility to get the bayes factor? Any hint is highly appreciated! Thank you very much.

    Best,
    Antonius


    This is my model:

    model {
    
      # choice Probabilities
      for (i in 1:ns){ #subject (total = 37)
        for (ii in 1:nDepri){ #Deprivation (total = 2)
          for (iii in 1:nTask){ #Task (total = 2)
            for (j in 1:nt[ii,iii,i]){ #Trial (total ~ 108 per Subject / Deprivation / Task)
    
              choice[ii,iii,i,j] ~ dbern(p[ii,iii,i,j,mdlIdx])
    
              # Model 1
              p[ii,iii,i,j,1] <- exp(sv[ii,iii,i,j]/(temp[ii,iii,i,j])) / ((exp(amount_now[ii,iii,i,j]/(temp[ii,iii,i,j]))) +                                                                         exp(sv[ii,iii,i,j]/(temp[ii,iii,i,j])))
              temp[ii,iii,i,j] <- (temp.slope[ii,iii,i]*trial[ii,iii,i,j]) + temp.intercept[ii,iii,i]
    
              # Model 2
              p[ii,iii,i,j,2] <- exp(sv[ii,iii,i,j]/(temp0[ii,iii,i,j])) / ((exp(amount_now[ii,iii,i,j]/(temp0[ii,iii,i,j]))) +                                                                         exp(sv[ii,iii,i,j]/(temp0[ii,iii,i,j])))
              temp0[ii,iii,i,j] <- (temp0.slope[i]*trial[ii,iii,i,j]) + temp0.intercept[i]
    
              #shared between both models
              sv[ii,iii,i,j] <- amount[ii,iii,i,j] / (1 + ((k[ii,iii,i])* delay_diff[ii,iii,i,j]))
    
            }
          }
        }
      }
    
      # Parameter for Subject 
      for (i in 1:ns){
        for (ii in 1:nDepri){
          for (iii in 1:nTask){
            k[ii,iii,i] ~ dnorm(k_mu[ii,iii],k_lambda[ii,iii])T(0.0001,)
            temp.slope[ii,iii,i] ~ dnorm(temp.slope_mu[ii,iii],temp.slope_lambda[ii,iii])
            temp.intercept[ii,iii,i] ~ dnorm(temp.intercept_mu[ii,iii],temp.intercept_lambda[ii,iii])T(.0001,)
    
          }
        }
    
        temp0.slope[i] ~ dnorm(temp0.slope_mu,temp0.slope_lambda)
        temp0.intercept[i] ~ dnorm(temp0.intercept_mu,temp0.intercept_lambda)T(0.0001,)
    
      }
    
      # Prior Gruppenverteilung BOTH MODELS
    
      k.m  <- 0.044
      k.sd  <- 0.0321
      k.prec <- 1/(k.sd^2)
    
      #1 non deprived - immediate
    
      k_mu[1,1] ~ dnorm(k.m,k.prec)T(0.0001,)
      k_lambda[1,1] ~ dgamma(.001,.001)T(.0001,)
    
      #2 non deprived - delayed
    
      k_mu[1,2] ~ dnorm(k.m,k.prec)T(0.0001,)
      k_lambda[1,2] ~ dgamma(.001,.001)T(.0001,)
    
      #3  deprived - immediate
    
      k_mu[2,1] ~ dnorm(k.m,k.prec)T(0.0001,)
      k_lambda[2,1] ~ dgamma(.001,.001)T(.0001,)
    
      #4  deprived - delayed
    
      k_mu[2,2] ~ dnorm(k.m,k.prec)T(0.0001,)
      k_lambda[2,2] ~ dgamma(.001,.001)T(.0001,)
    
      # Prior Gruppenverteilung - Model 1
    
      #1 non deprived - immediate
    
      temp.slope_mu[1,1] ~ dnorm(temp.slope.m[1],temp.slope.prec[1])
      temp.slope_lambda[1,1] ~ dgamma(.001,.001)T(.001,)
    
      temp.intercept_mu[1,1] ~ dnorm(temp.intercept.m[1],temp.intercept.prec[1])T(.0001,)
      temp.intercept_lambda[1,1] ~ dgamma(.001,.001)T(.001,)
    
      #2 non deprived - delayed
    
      temp.slope_mu[1,2] ~ dnorm(temp.slope.m[2],temp.slope.prec[2])
      temp.slope_lambda[1,2] ~ dgamma(.001,.001)T(.001,)
    
      temp.intercept_mu[1,2] ~ dnorm(temp.intercept.m[2],temp.intercept.prec[2])T(.0001,)
      temp.intercept_lambda[1,2] ~ dgamma(.001,.001)T(.001,)
    
      #3  deprived - immediate
    
      temp.slope_mu[2,1] ~ dnorm(temp.slope.m[3],temp.slope.prec[3])
      temp.slope_lambda[2,1] ~ dgamma(.001,.001)T(.001,)
    
      temp.intercept_mu[2,1] ~ dnorm(temp.intercept.m[3],temp.intercept.prec[3])T(.0001,)
      temp.intercept_lambda[2,1] ~ dgamma(.001,.001)T(.001,)
    
      #4  deprived - delayed
    
      temp.slope_mu[2,2] ~ dnorm(temp.slope.m[4],temp.slope.prec[4])
      temp.slope_lambda[2,2] ~ dgamma(.001,.001)T(.001,)
    
      temp.intercept_mu[2,2] ~ dnorm(temp.intercept.m[4],temp.intercept.prec[4])T(.0001,)
      temp.intercept_lambda[2,2] ~ dgamma(.001,.001)T(.001,)
    
      #acutal priors
    
      temp.slope.m[1]  <- 0
      temp.slope.m[2]  <- 0
      temp.slope.m[3]  <- 0
      temp.slope.m[4]  <- 0
    
      temp.slope.prec[1] <- 1/(0.01^2)
      temp.slope.prec[2] <- 1/(0.01^2)
      temp.slope.prec[3] <- 1/(0.01^2)
      temp.slope.prec[4] <- 1/(0.01^2)
    
      #model 1
      temp.intercept.m[1]  <- 3
      temp.intercept.m[2]  <- 3
      temp.intercept.m[3]  <- 3
      temp.intercept.m[4]  <- 3
    
      temp.intercept.prec[1] <- 1/(0.1^2)
      temp.intercept.prec[2] <- 1/(0.1^2)
      temp.intercept.prec[3] <- 1/(0.1^2)
      temp.intercept.prec[4] <- 1/(0.1^2)
    
      # Prior group distrebution - Model 2
    
      #1 non-deprived
      temp0.slope_mu ~ dnorm(temp0.slope.m,temp0.slope.prec)
      temp0.slope_lambda ~ dgamma(.001,.001)T(.0001,)
      temp0.slope_sigma <- 1/sqrt(temp0.slope_lambda)
    
      temp0.intercept_mu ~ dnorm(temp0.intercept.m,temp0.intercept.prec)T(.0001,)
      temp0.intercept_lambda ~ dgamma(.001,.001)T(.0001,)
      temp0.intercept_sigma <- 1/sqrt(temp0.intercept_lambda)
    
      temp0.slope.m  <- 0
      temp0.slope.prec <- 1/(0.01^2)
      temp0.intercept.m  <- 3
      temp0.intercept.prec <- 1/(0.1^2)
    
      # Prior on Model Index
      mdlIdx ~ dcat( modelProb[] ) # Index 1 = model1 | index 2 = model2
    
      # Hyperprior on model index:
      modelProb[1] <- prior1
      modelProb[2] <- 1-prior1
    
    }
    
     
  • Martyn Plummer
    Martyn Plummer
    2013-10-10

    To get jumping between models you need to use the approach of Carlin and Chib: set up pseudo-priors that approximate the posterior distribution of each model when the model is not active. This requires you to do some pilot runs for each model so you can approximate the posterior.

    A worked example is given here:

    http://www.mrc-bsu.cam.ac.uk/bugs/winbugs/examples/pines.txt

    The reason it won't jump without the pseudo-priors is that when the model is not active it's parameters are sampled from their prior distribution and are typically very far from the set of values that are supported by the data. Hence a jump that makes the inactive model active is highly unlikely.

     
  • Hi Martyn,
    thank you for your quick and well-informed reply! I'll try to implement the pseudo prior approach now. Thank you!

     
  • Hi,
    I set up a model using pseudo priors on single subject level. Unfortunately, even with a model 1 prior <0.001 there is not a single jump to the second model. For monitoring reasons, I tracked the deviance, too. It attracted my attention that the deviance chains are not mixing very well in 60000 iterations. Is there anything, that I could do optimized this model and sample the model index or calculate the bayes factor?

    Any advice would be very helpful and highly appreciated!

    Thanks a lot!

    Best,
    Antonius