#Example:hierachical logistic regression model for rat tumor rates 
#Gelman and Rubin book, pag 291
#
fun_function(x, y)
{(x -y+5)^2+5}
 test_nlminb(c(0,0),fun)

make.data_ function() {
  tmp _ cbind(c(4,4,2),c(14,34,34),c(0,1,2)) 
  tarone _ matrix(scan("tarone.txt"),byrow=T,ncol=2)
  yy0    _ tarone[,1]
  nn0   _ tarone[,2]
 yy _ tmp[,1]       #number of tumors
  nn _ tmp[,2] #total number of births
  xx _ tmp[,3]       #dose 
  sf _ cbind(yy,nn-yy)
  sf0 _ cbind(sum(yy0),sum(nn0)-sum(yy0)) #historical data
  sf0_rbind(sf0,sf[2,],sf[3,])
#crude estimate of tau
for( i in 1 : length(yy0)){
if (yy0[i] == 0) yy0[i] _ .5}
tau.hat_sqrt(var(log( (yy0/nn0)/(1-yy0/nn0)))) 
 return(tmp,sf,xx,sf0,tau.hat)
}

separate_function(){
  DD_make.data() 
  glm.out        _ glm(sf ~ xx, family=binomial, data=DD)
  return(glm.out)
}

pooling_function(){
  DD_make.data() 
  glm.out        _ glm(sf0 ~ xx, family=binomial, data=DD)
  return(glm.out)
}

##non hierarchical logistic regression
glm.separate_separate()
glm.pooling_pooling()


#find conditional mode

logpost_function(gamma,tau=1){
  yy_DD$tmp[,1]
  nn_DD$tmp[,2]
  xx_DD$xx 
  tarone _ matrix(scan("tarone.txt"),byrow=T,ncol=2)
  yy0   _ tarone[,1]
  nn0   _ tarone[,2]
  J_length(yy0)-1
  alpha_gamma[(1:(J+1))]
  beta_gamma[J+2]
  mu_gamma[J+3]
   - (alpha[1:J]%*%yy0[1:J]-nn0[1:J]%*%log(1+exp(alpha[1:J]))+(J+1)*log(tau)
	  - (1/(2*tau^2))*sum(alpha[1:J]-mu)^2
	  + alpha[J+1]*sum(yy)
	  + beta*xx%*%yy
	  - nn%*%log(1+exp(alpha[J+1]+beta*xx)))
          #min(-logposterior)=max(logposterio)
}


    mode_ nlminb(start=c(rep(0,71),0,0),obj=logpost)