library(rjags)

dataimp = read.delim("http://web.as.uky.edu/statistics/users/pbreheny/701/data/brunner.txt")
data = list(p.model=c(0.2,0.2,0.2,0.2,0.2),
            returned=dataimp$returned,
            S=dataimp$S, # centered
            L=dataimp$L) # standardized

SL = data$L*(data$S>-1) # data$S < -1 are 0's
SLpar = c(mean(SL),sd(SL))
data$SL = (SL-SLpar[1])/SLpar[2]

data$n = length(data$returned)
data$sexpred = -1.41494 
data$lenpred = 1.5
data$sexlenpred = ((data$sexpred>-1)*data$lenpred-SLpar[1])/SLpar[2]

nchain = 3; nadapt = 5000; niter = 10000
m1 = jags.model("http://web.as.uky.edu/statistics/users/pbreheny/701/S13/notes/2-21.txt", data = data, n.chains = nchain, n.adapt = nadapt)
out1 = jags.samples(m1,c("beta0","beta1","beta2","beta12","mod","pred.prob"), n.iter = niter)

#### First -- find model probs
modprob = table(out1$mod)/(niter*nchain)

#### Next -- find BF
nmod = 5
bf = matrix(NA,nmod,nmod)
for(j in 1:nmod){
  bf[j,] = modprob[j]/modprob
}

#### find predictions (overall) and within each model
whchmod = vector('list',nmod)
densmod = vector('list',nmod)
allx = ally = NULL
modvec = c(out1$mod); predvec = c(out1$pred.prob)
for(j in 1:nmod){
  whchmod[[j]] = which(modvec==j)
  densmod[[j]] = density(predvec[whchmod[[j]]])
  allx = c(allx,densmod[[j]]$x)
  ally = c(ally,densmod[[j]]$y)
}
plot(NULL,NULL,xlim=range(allx),ylim=range(ally),xlab="Predicted probability of return",yaxt="n",ylab="") #xlim=range(allx)
for(j in 1:nmod){
  lines(densmod[[j]]$x,densmod[[j]]$y,col=j+1)
}
denall = density(out1$pred.prob)
lines(denall$x,denall$y,lwd=2)
legend("topright",c("Model 1", "Model 2", "Model 3","Model 4", "Model 5", "Model Averaged"),col=c(2:6,1),lwd=c(rep(1,5),2),lty=1)

#### Get traceplot of beta2 -- not perfect, but good enough
ntot = niter*nchain
coluse = modvec>=3
lwduse = (modvec>=3)*2
plot(NULL,NULL,xlim=c(0,ntot),ylim=range(out1$beta2),xlab="Iteration",ylab="Value")
bet2vec = c(out1$beta2)
lines(c(1,1.5),c(bet2vec[1],mean(bet2vec[1:2])),col=coluse[1]+1,lwd=lwduse[1]+1)
for(i in 2:(ntot-1)){
  lines(c(i-0.5,i,i+0.5),c(mean(bet2vec[(i-1):i]),bet2vec[i],mean(bet2vec[i:(i+1)])),col=coluse[i]+1,lwd=lwduse[i]+1)
}
lines(c(ntot-1,ntot),c(mean(bet2vec[(ntot-1):(ntot)]),bet2vec[ntot]),col=coluse[ntot]+1,lwd=lwduse[ntot]+1)