#################################################
# R-code used to call WinBUGS
# ex: Logistic regression: Dose-response study
# SIMULATION STUDY
#################################################


#------------------------------------------------
# Call Package
#------------------------------------------------
 library(R2WinBUGS) # calls WinBUGS from R

#------------------------------------------------
# Fixed Parameter values
#------------------------------------------------
  N.sim=5
  k=6
  n=rep(50,k)
  dose=c(1,15,25,30,40,50)
  beta.true=c(-10, 0.4)
  prob.true=1/(1+exp(-beta.true[1]-beta.true[2]*dose))
  
#------------------------------------------------
#  Define Simulation Result Matrices
#------------------------------------------------
  p=2
  results_MEAN=matrix(0,nrow=N.sim,ncol=p)
  results_MEDIAN=matrix(0,nrow=N.sim,ncol=p)
  results_SD=matrix(0,nrow=N.sim,ncol=p)
  results_QU2.5=matrix(0,nrow=N.sim,ncol=p)
  results_QU97.5=matrix(0,nrow=N.sim,ncol=p)

############################################################
################# Start Simulation #########################
############################################################


start.time= proc.time()     # to measure run time for simulation

for(i in 1:N.sim){

#------------------------------------------------
# Generate Data
#------------------------------------------------

 x=rbinom(k,size=n,prob=prob.true)

#------------------------------------------------
# Inputs to the bugs function
#------------------------------------------------
  BUGScode<-file.path("E:/Codes/R-WinBUGS/logitmodel.txt")
##THE ABOVE FILE NAME PATH SHOULD BE CHANGED BUT FILENAME SHOULD BE THE SAME

    # usually you want more runs than this!!!
    # (but this will run faster during class)
  runs<-1100     # total number of runs (including burnins)
  burn<-100     # number of burnins
  nthin<-1
  nchains<-3
  debug<-TRUE
  data<-list("n", "x", "dose", "k")
  keepers<-list("beta")
  geninits1<-list(beta=c(0,0))
  geninits2<-list(beta=c(-1,-0.5))
  geninits3<-list(beta=c(1,0.5))

#------------------------------------------------
# Call Winbugs
#------------------------------------------------

  fit<-bugs(
    model.file=BUGScode,
    data=data,
    inits = list(geninits1,geninits2,geninits3),
    parameters.to.save = keepers,
    n.chains=nchains,
    n.iter=runs,
    n.burnin=burn,
    n.thin=nthin,
    #codaPkg=TRUE,  #output diagnostic information
    bugs.directory="E:/WinBUGS14/",
#THE ABOVE FILE NAME PATH SHOULD BE CHANGED (typically it is C:/Program Files/WinBUGS14)

    #debug=debug,
    clearWD=TRUE,
    working.directory=("E:/Codes/R-WinBUGS/")
   #THE ABOVE FILE NAME PATH SHOULD BE CHANGED

  )

#------------------------------------------------
# Place Summary Statistics in Matrices
#------------------------------------------------
  results_MEAN[i,1:p]=  fit$summary[1:p,1]
  results_SD[i,1:p] =   fit$summary[1:p,2]
  results_QU2.5[i,1:p]= fit$summary[1:p,3]
  results_MEDIAN[i,1:p]=fit$summary[1:p,5]
  results_QU97.5[i,1:p]=fit$summary[1:p,7]

  #print(i)
  }

############################################################
################### End Simulation #########################
############################################################

#------------------------------------------------
# Measure run time for simulation
#------------------------------------------------

  run_time=proc.time()-start.time #measured in seconds
  run_minutes=run_time[3]/60
  run_minutes

#------------------------------------------------
# Compute Summary Statistics
#------------------------------------------------

  MEAN = apply(results_MEDIAN,2,mean)   # mean
  TRUTH=beta.true
  BIAS = MEAN-TRUTH                 # bias
  MCSE = apply(results_MEDIAN,2,sd)     # monte carlo se (sd(median))
  ESE  = apply(results_SD,2,mean)       # est. se (mean(sd))
  MSE=MEAN^2+MCSE^2
  COV=c()
  for(i in 1:p){
COV[i]=sum(as.numeric(results_QU2.5[,i]<=TRUTH[i] & results_QU97.5[,i] >=TRUTH[i]))/N.sim}


  summary=t(rbind(TRUTH, MEAN,BIAS,MCSE,ESE,MSE,COV))
  rownames(summary)<-c("beta1","beta2")
  colnames(summary)<-c("TRUE", "MEAN","BIAS","MCSE","ESE","MSE","COV")

  SUM=round(summary,3)

  print(SUM)

#------------------------------------------------
#  Boxplots
#------------------------------------------------

  boxplot(results_MEDIAN[,1], results_MEDIAN[,2],
        names=c("beta1","beta2"),main=paste("#MC runs=",N.sim))

  abline(h=TRUTH[1],lty=2,col="red")
  abline(h=TRUTH[2],lty=2,col="red")
  
  cat(c("Elasped time in minutes:",round(run_minutes,2)),fill=T)