#Read in counties.txt
#R uses forward slashes for directory structures rather than backward slashes

county<-read.table('c:/r/counties.txt',header=TRUE)


#Generate MLE of county specific rates and exact confidence intervals
mle.exact<-pois.exact(county$cases,county$pop)
#Results are saved in mle.exact, take a look
mle.exact


#Run MCMC using openbugs
#define the data by generating vectors for cases and population
cases<-county$cases
pop<-county$pop
#This step combines the data into a format that BUGS can use
data<-list("cases","pop")
#Now define the initial values...in this case theta needs initial values
#there are 100 values of theta and we set them all to 0 initially
#rep(0,100) generates a vector of 100 zeros
inits <- function() {list (theta=rep(0,100)) }
#Now we tell bugs which parameters we want to keep track of (theta, and et=exp(theta))
parameters<-c("theta","lambda")
#Now, run bugs passing in the data, inits, parameters and the bugs program counties.bug.  run 1 chain for 5000 iterations,discarding the first 1000.  We're using openbugs, but you could use winbugs
county.sim<-bugs(data,inits,parameters,"counties.bug",n.chain=1,n.iter=5000,n.burnin=1000,n.thin=1,program="openbugs")
#this step attaches the results from bugs so you can look at them
attach.bugs (county.sim)
#Look at the results
print(county.sim)
#look at convergence (this is for the first theta)
plot(theta[,1])


#Generate the Plots
#get the order of coefficient magnitudes so we can sort by that
oeff<-order(mle.exact$rate,decreasing=TRUE)
plot(county[,3],mle.exact$upper[oeff],ylim=range(c(mle.exact$lower,mle.exact$upper)),col="black",type='n',xlab="County", ylab="Breast Cancer Rate/ 100,000")
segments(county[,3], mle.exact$lower[oeff], county[,3], mle.exact$upper[oeff]) 
points(county[,3],mle.exact$rate[oeff],col="black",pch=19)
mcrude<-mean(mle.exact$rate)
lines(c(0,100),c(mcrude,mcrude),lty=1,col="green")