Stratified_NevadaThyroiditis<-function(fractionBerk)

##########################################################################################################
##########################################################################################################
#
#Author: Ciprian Craniceanu
#e-mail: ccrainic@jhsph.edu
#Last revised: February 10, 2008 by Carroll
#
#
##########################################################################################################
##########################################################################################################
#
#Description: 
#	1. Fits an excess relative risk logistic model with measurement error, using the Nevada Test Site Thyroid Disease Study Data. 
#     2. WinBUGS model is stored in file CMlogisticNevadaingham.txt
#	3. The mean and standard deviations of the outcome model parameters are stored and returned 
#	4. You can manipulate how much of the measurement error is Berkson, and how much is classical. You do this
#		via the input argument fractionBerk, which is the fraction (between 0 and 1) of the
#		measurement error that is Berkson
#
##########################################################################################################
##########################################################################################################
#
#Input:
#
#
##########################################################################################################
##########################################################################################################
#
#Output:
#	1. Given in the structure called results at the end of the program
#
##########################################################################################################
##########################################################################################################
#
# After you source this program (source("Stratified_NevadaThyroiditis.R")) the Call as a follows:
#
#                             Stratified_NevadaThyroiditis(fractionBerk)
#
#
##########################################################################################################
##########################################################################################################

{#begin function

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

#Define the constants of the data and model
Nobservations=2491 	#Total number of observations. This has to be exactly correct!
#Nobservations= 500 	#Total number of observations. This has to be exactly correct!

#Load and assign data
#data.Nevada<-as.matrix(read.table(file="C:\\pctex_Raymond_Carroll_directory\\papers\\utah_2005_directory\\final_pha#se3_data\\fake_utah_stratified.txt"))
data.Nevada<-as.matrix(read.table(file="enar_utah.txt"))

Age     <- data.Nevada[,2]
Sex     <- data.Nevada[,3]
Sex     <- Sex - 1            	# The women have value = 1
geomean <- data.Nevada[,4]    	# The geometric mean
geostd  <- data.Nevada[,5]    	# The geomatrix standard deviation
Y       <- data.Nevada[,1]   		# The binary response
W       <- log(0.01 + geomean)      # The values of W used
stratum <- data.Nevada[,6]	      # The stratum

sigmasqtotal <- log(geostd) ^ 2;					# The total uncertainty
#fractionBerk ; 								# The fraction of the uncertainty that is Berkson
sigmasqB     <- sigmasqtotal * fractionBerk     		# The amount of Berkson uncertainty
sigmasqC     <- sigmasqtotal * (1- fractionBerk)       	# The amount of classical uncertainty
tauB       <- 1 / sigmasqB                                  # WinBUGS is strange and insists on using the inverse 
                                                            # of variances
tauC       <- 1 / sigmasqC                                  # Double Strange

muL        <- mean(W)							# The mean of W as a starting value for the mean 
                                                            # of the latent intermediate L
sigmasqL   <- var(W) - mean(sigmasqC)				# The starting value of the latent intermediate L
tauL       <- 1 / sigmasqL

# There is a bit of a tricky issue in getting starting values for the
# latent intermediate value and also the latent X. I use regression
# calibration, i.e., I generated L from the conditional distribution of L given W,
# and then I generated X from the Berkson part of the model.
#
# It is possible to generate the initial L from the normal conditional
# distribution of L given W.
attenuation  <- sigmasqL / var(W)
inits.L      <- muL + (attenuation * (W-muL))
inits.L      <- rnorm(Nobservations,mean=inits.L,sd = sqrt(attenuation * sigmasqC))
inits.X      <- rnorm(Nobservations, mean=inits.L, sd=sqrt(sigmasqB)) 


##################################################################################################
##################################################################################################
#
#
# Define data, initial values and parameters to be monitored in the Bayesian analysis 
# WinBUGS insists that anything passed is also used
#
##################################################################################################
##################################################################################################

#Data
	data<-list("Y","W","Sex","Nobservations","sigmasqB","tauB","tauC","Age","stratum")



#Initial values
# beta has four entries
#      intercept
#      female effect
#      Age effect
#      Excess relative risk
# alpha has the mean of the latent intermediate for each stratum
# gamma has the inverse of the variance for each stratum
# 
#     Variance of the latent intermediate
	inits<-function(){
	list(
	   beta = c(-4.50,1.75,0.53,6.00),
        alpha = c(muL,muL,muL,muL,muL),
        gamma = c(1/tauL,1/tauL,1/tauL,1/tauL,1/tauL),
	      X = inits.X,
            L = inits.L)
	}

#Parameters
parameters<-list("beta","alpha","gamma")



# Run the Bayesian analysis of the Nevadaingham data. 
# You determine the burn in
# You determine the number of simulations from the posterior distribution of the parameters given the data.
# You determine through nthin to retain every nthin simulated sample

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

results<-fit$sims.matrix
beta  <- results[,1:4]
alpha <- results[,5:9]
gamma <- results[,10:14]


plot(results[,4],type="l")

output.parameters<-cbind(alpha,beta,gamma)

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

#########################################################################
# Now write out your results for the various parameters.
# It is useful to name the file something that can be remebered :)
#########################################################################

write.table(results,file="Raystable_test01_fracberk01.txt")


}#end function