#################################################################################
# Bayesian model selection and statistical modeling 
# Chapman & Hall/CRC Taylor and Francis Group 
#
# Chapter7: Section7.2.1
#
# Hierarchical Bayesian modeling for logistic regression
#
# Author : Tomohiro Ando
#
################################################################################
library(MCMCpack)

# Dutch voting behavior data

D <- matrix(scan("Defaultdata.txt"),ncol=4,byrow=T)
D <- as.data.frame(D)
names(D) <- c("ROA","logSales","CF-Sales ratio","Ratings")
n <- nrow(D)


library(bayesm)

R <- 6000

nvar=2    # number of coefficients
nlgt=50   # number of cross-sectional units
nobs=5    # number of observations per unit
nz=3      # number of regressors in mixing distribution

## set hyper-parameters
##     B=ZDelta + U  

Z=matrix(c(rep(1,nlgt),runif(nlgt,min=20,max=50),trunc(runif(nlgt,min=0.0001,max=1.9999))),nrow=nlgt,ncol=nz)
Delta=matrix(c(-2,-1,0,-1,1,-.5),nrow=nz,ncol=nvar)
iota=matrix(1,nrow=nvar,ncol=1)
Vbeta=diag(nvar)+0.5*iota%*%t(iota)

## simulate data
lgtdata=NULL

for (i in 1:nlgt) 
{ beta=t(Delta)%*%Z[i,]+as.vector(t(chol(Vbeta))%*%rnorm(nvar))
  X=matrix(runif(nobs*nvar),nrow=nobs,ncol=nvar)
  prob=exp(X%*%beta)/(1+exp(X%*%beta)) 
  unif=runif(nobs,0,1)
  y=ifelse(unif<prob,1,0)
  lgtdata[[i]]=list(y=y,X=X,beta=beta)
}

out=rhierBinLogit(
Data=list(lgtdata=lgtdata,Z=Z),
Prior=list(nu=10),
Mcmc=list(R=R,keep=5))


PD <- max((out$llike)[-(1:100)])-mean((out$llike)[-(1:100)])
DIC <- -2*mean(out$llike)+2*PD

print(DIC)