dat<-read.csv("C:\\Users\\bjreich.ST\\Documents\\ST740\\Votes08.csv")
dat.adj<-read.csv("C:\\Users\\bjreich.ST\\Documents\\ST740\\NCADJ.csv")

county.names<-dat.adj[,1]
ADJ<-dat.adj[,-1]
dimnames(ADJ)<-list(county.names,county.names)

library(maps)

plot.nc<-function(data,title=""){
   d<-data[c(1:26,27,27,27:100)] #county 27 has three polygons!
   map("county","north carolina",fill=T,col=d)
   title(title)
}
   
logit<-function(x){exp(x)/(1+exp(x))}

Logit.CAR<-function(n,y,ADJ,n.samples=5000,burn=1000){
    #fit the model:
    #  y~binom(n,logit(theta))
    #  theta~CAR

    n.s<-nrow(ADJ)
    m<-apply(ADJ,2,sum)
    Q<-as.matrix(diag(m)-ADJ)
    ADJ<-ADJ==1

    #initial values:
    theta<-rep(0,n.s)
    tau<-1

    acc<-keeptau<-rep(0,n.samples)
    keeptheta<-matrix(0,n.samples,n.s)

    for(i in 1:n.samples){
      #update theta:
      for(s in 1:n.s){
        cantheta<-theta
        cantheta[s]<-rnorm(1,theta[s],.5)
        R<-dbinom(y[s],n[s],logit(cantheta[s]),log=T)-
           dbinom(y[s],n[s],logit(theta[s]),log=T)+
           dnorm(cantheta[s],mean(theta[ADJ[s,]]),1/sqrt(m[s]*tau),log=T)-
           dnorm(theta[s],mean(theta[ADJ[s,]]),1/sqrt(m[s]*tau),log=T)
        if(runif(1)<exp(R)){
            theta<-cantheta
            acc[i]<-acc[i]+1/n.s
        }
      }
      tau<-rgamma(1,(n.s-1)/2+0.01,rate=t(theta)%*%Q%*%theta/2+0.01)

      keeptheta[i,]<-theta
      keeptau[i]<-tau
      if(i%%50==0){plot.nc(ifelse(theta>0,4,2),i)}
    }

p<-logit(keeptheta[burn:n.samples,])

list(p.samples=p,p.mn=apply(p,2,mean),p.sd=apply(p,2,sd),
     MH.acc=acc,tau=keeptau)}

#The full data set:
y.total<-dat$Obama
n.total<-dat$Total
p.total<-y.total/n.total

#Generate a random subset:
n<-rbinom(100,n.total,0.001)
y<-rbinom(100,n,p.total)

#Fit the CAR model:
fit<-Logit.CAR(y=y,n=n,ADJ=ADJ)

#plot results
plot.nc(ifelse(p.total>0.5,4,2),"All data")
plot.nc(ifelse((y+1)/(n+1)>0.5,4,2),"Subset of  data")
plot.nc(ifelse(fit$p.mn>0.5,4,2),"All data")
hist(fit$p.sd)

#Compare the fits of the CAR model and raw data estimate:
MAD.spatial<-100*mean(abs(p.total-fit$p.mn))
MAD.raw<-100*mean(abs(p.total-(y+1)/(n+1)))
print("MAD, spatial v non-spatial")
print(MAD.spatial)
print(MAD.raw)