MainprogramlogisticFramingham<-function()

##########################################################################################################
##########################################################################################################
#
#Author: Ciprian Craniceanu
#e-mail: ccrainic@jhsph.edu
#Last revised: December 29, 2005
#
#
##########################################################################################################
##########################################################################################################
#
#Description: 
#	1. Fits a logistic linear model with measurement error. WinBUGS model is stored in file 
#	   logisticFramingham.txt
#	3. The mean and standard deviations of the outcome model parameters are stored and returned 
#	4. Results are contrasted with the naive analysis where the mean of W's is used as the true X
#
##########################################################################################################
##########################################################################################################
#
#Input:
#
#
##########################################################################################################
##########################################################################################################
#
#Output:
#	1. output.Bayes:	An Nsimulations*6 matrix containing the posterior mean and standard 
#				deviations of the beta parameters for the Bayes analysis. 
#
##########################################################################################################
##########################################################################################################
#
#Call:
#
#MainprogramlogisticFramingham()
#
#
##########################################################################################################
##########################################################################################################

{#begin function

#Uses R2WinBUGS as an interface between R and WinBUGS
library(R2WinBUGS)

#Define the constants of the data and model
Nobservations=641 		#Total number of observations 
Nreplications=2			#Total number of replications
nalphas=2			#Number of alpha (exposure model) parameters
nbetas=3			#Number of beta (outcome model) parameters

#Load and assign data
data.Fram<-as.matrix(read.table(file="Framinghamdata.txt"))
W<-data.Fram[,1:2]
Z<-data.Fram[,3]
Y<-data.Fram[,4]
inits.X<-(W[,1]+W[,2])/4



output.Bayes<-rep(0,6)
dim(output.Bayes)<-c(1,6)


##################################################################################################
##################################################################################################
#
#
#Define data, initial values and parameters to be monitored in the Bayesian analysis 
#
##################################################################################################
##################################################################################################

#Data
	data<-list("Y","W","Z","nalphas","nbetas","Nobservations","Nreplications")



#Initial values
	inits<-function(){
	list(
	alpha=c(4.42,-0.02),beta=c(-10.10,1.76,0.38),
	taux=20,lambda=0.3,X=inits.X)
	}

#Parameters
parameters<-list("alpha","beta","sigmax2","sigmau2","lambda")

#Run the Bayesian analysis of the Framingham data. Use 10000 burn-in and 300000
#simulations from the posterior distribution of the parameters given the data.
#Retain only every 60th simulated sample

date()
fit<- bugs(data, inits, parameters, model.file = "logisticFramingham.txt",
         n.chains = 1, n.iter = 1000, n.burnin = 10,
         n.thin = 10,debug = FALSE, DIC = FALSE, digits = 5, codaPkg = FALSE,
         bugs.directory = "c:/Program Files/WinBUGS14/")
date()
attach.all(fit)

output.Bayes[seq(1,6,by=2)]=fit$mean$beta
output.Bayes[seq(2,6,by=2)]<-fit$sd$beta


write.table(output.Bayes,file="Bayesbetaparametersoutput.txt")
return(list(output.Bayes))

output.parameters<-cbind(alpha,beta,sigmax2,sigmau2,lambda)

write.table(output.parameters,file="Bayeshistory.txt")

}#end function