## transdimensional model switch index is not jumping at all document.SUBSCRIPTION_OPTIONS = { "thing": "thread", "subscribed": false, "url": "subscribe", "icon": { "css": "fa fa-envelope-o" } };

Help
2013-10-09
2013-10-30
• Antonius Wiehler - 2013-10-09

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 (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){
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 - 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.

• Antonius Wiehler - 2013-10-15

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

• Antonius Wiehler - 2013-10-30

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?

Thanks a lot!

Best,
Antonius