0

I'm incorporating a Bayesian Model Averageing(BMA) approach in my research and strapped in trapped in the estimated of Pr(theta|D). Professor John K. Kruschke(2014)'s book in chapter 10 offers an example to compute model posterior probability using hierarchical MCMC, but I do not know how to set the index value of model(m) using this method to linear regression. Here is the code

library(BMS)
data(attitude)
#show the model posterior probability of four combination of variable compaints and #privileges
bma4c = bms(rating ~complaints + privileges , data = attitude)
pmp.bma(bma4c)

I want use hierarchical MCMC method to get the similar results.

# The date
library(rjags)
dataList <- list( N = nrow(attitude),
                  y =attitude$rating,
                  x = attitude$complaints,
                  z = attitude$privileges)


#The model

modelString = "
model {
for ( i in 1:N ) {
y[i] ~ dnorm( mu[i,m] , tau )
mu[i,1] <- beta0 + beta1 *x 
mu[i,2] <- beta0  
mu[i,3] <- beta0 + beta2 *z 
mu[i,4] <- beta0 + beta2 *z +beta1 *x 

}
tau <- 1/100

mu[] <- equals(m,1)*mu[,1] + equals(m,2)*mu[,2] + equals(m,3)*mu[,3] + equals(m,4)*mu[,4]
# this sentence cause the error  <Cannot evaluate subset expression for mu>

beta0 ~ dnorm(0 , 1.0E-12)
beta1 ~ dnorm(0 , 1.0E-12)
beta2 ~ dnorm(0 , 1.0E-12)

m ~ dcat(mPrior[])
mPrior[1] <- 0.25
mPrior[2] <- 0.25
mPrior[3] <- 0.25
mPrior[4] <- 0.25
}" 


 writeLines( modelString , con="TEMPmodel.txt" )

# RUN THE CHAINS.

parameters = c("theta","m") 
adaptSteps = 1000             # Number of steps to "tune" the samplers.
burnInSteps = 1000           # Number of steps to "burn-in" the samplers.
nChains = 4                   # Number of chains to run.
numSavedSteps=50000          # Total number of steps in chains to save.
thinSteps=1                   # Number of steps to "thin" (1=keep every step).
nPerChain = ceiling( ( numSavedSteps * thinSteps ) / nChains ) # Steps per chain.
# Create, initialize, and adapt the model:
jagsModel = jags.model( "TEMPmodel.txt" , data=dataList , 
                        n.chains=nChains , n.adapt=adaptSteps )
# Burn-in:
update( jagsModel , n.iter=burnInSteps )
# The saved MCMC chain:
codaSamples = coda.samples( jagsModel , variable.names=parameters , 
                            n.iter=nPerChain , thin=thinSteps )
mcmcMat = as.matrix( codaSamples , chains=TRUE )
m = mcmcMat[,"m"]
# compute the model posterior probability
pM1 = sum( m == 1 ) / length( m )
pM2 =  sum( m == 2 ) / length( m )
pM3 =  sum( m == 3 ) / length( m )
pM4 =  sum( m == 4 ) / length( m )

Does anyone know how to modify this code to get the relative model probability?

ilanman
  • 4,787
Jieyu He
  • 3

2 Answers2

0

Two different things to think about: How to code the models in a logically convenient way, and how to get MCMC to sample efficiently from the models. Even if you solve the first, it might be difficult to solve the second.

First, it appears that your model is just multiple regression with two predictor variables that you call x and z. Your different models consider including or excluding x and z. This situation is typically called "variable selection," and this exact situation is addressed extensively in Section 18.4 of the 2nd edition of Doing Bayesian Data Analysis. The idea is that every predictor variables is provided with a multiplicative inclusion parameter delta[i] that can take on values of 0 or 1. A version of the analysis was discussed on my blog before the 2nd Edition was published, here: http://doingbayesiandataanalysis.blogspot.com/2014/01/bayesian-variable-selection-in-multiple.html A follow-up regarding $R^2$ was posted here: http://doingbayesiandataanalysis.blogspot.com/2016/07/bayesian-variable-selection-in-multiple.html

Second, there's the issue of getting MCMC to efficiently sample from the four discrete models. The problem is that MCMC can get stuck and not jump between models well. The problem sometimes can be ameliorated by using pseudo-priors as discussed in DBDA2E. With a small set of models as in the present case, you can probably get nice mixing of the chains and a useful analysis. With larger sets of models, this method becomes unwieldy and ineffective.

John K. Kruschke
  • 2,361
  • Dear prof John, I written another code and got the relative model probability, the diagnosis of MCMC indicate the model index value jump between model properly, but such result did not observed if we treat the different combination of variable as different model, MCMC not jump between models(all model index value equal to 1) even I set the pseudo prior, why it not work? – Jieyu He Dec 29 '16 at 01:33
0

Thanks to Martyn Plummer and John K. Kruschke's guided, the MCMC can jump between model as fellow

model{
  for(i in 1:Ntotal){
    # log -logistic model
    theta[i,2] <- a2 + (1-a2)/(1+exp(-b2-c2*log(dose[i])))
    # log - probit  model
    probit(mid[i])  <- b1 + c1*log(dose[i])
    theta[i,1] <- a1 +(1-a1)*mid[i]

    mu[i]  <- theta[i,m]

    effect[i] ~ dbin(mu[i] ,N[i])
  }

  # Priors for model 2

  mu.a2 <- c(0.091, 0)
  mu.b2 <- c(-4.11, 0)
  mu.c2 <- c(1.61, 0)

  tau.a2 <- c(1957, 0.001)
  tau.b2 <- c(3.64, 0.001)
  tau.c2 <- c(22.79, 0.001)

  a2 ~ dnorm(mu.a2[m], tau.a2[m])
  b2 ~ dnorm(mu.b2[m], tau.b2[m])
  c2 ~ dnorm(mu.c2[m], tau.c2[m])

  # Priors for model 1

  mu.a1 <- c(0, 0.102)
  mu.b1 <- c(0, -2.52)
  mu.c1 <- c(0, 0.981)

  tau.a1 <- c(0.001, 2443)
  tau.b1 <- c(0.001, 11.21)
  tau.c1 <- c(0.001, 65.24)

  a1 ~ dnorm(mu.a1[m], tau.a1[m])
  b1 ~ dnorm(mu.b1[m], tau.b1[m])
  c1 ~ dnorm(mu.c1[m], tau.c1[m])

  ##index value
  m ~ dcat(mprior[])
  mprior[1] <-0.5
  mprior[2] <-0.5
}

But the calculated post-weight is not appropriate, more fine tuning and sampling method may need.

Jieyu He
  • 3