########################################
##  Analyze Data
##  josh clinton
##
##  princeton university
##  Clinton, Jackman and Rivers: "The Most Liberal Senator"
########################################

###############################################################################################
#  National Journal Votes Only
#  No Bush


idealout <- scan(file="C:/Documents and Settings/clinton/My Documents/ideal/Nat Journal/Key/src/x.out")
NJ.key <- matrix(idealout,ncol=101,byrow=T)
NJ.key <- NJ.key[-seq(1,200),-1]      #   drop start values + burn-in (200,000), iter count

xbar <- apply(NJ.key,2,function(x)quantile(x,c(.025,.50,.975)))
xbar <- t(xbar)
xbar[,2] <- apply(NJ.key,2,mean)


name <- scan(file="C:/Documents and Settings/clinton/My Documents/Projects/Nat. Journal/justname.txt",
                   what=list(name=""),
                   sep="\n")$name
                   
dimnames(xbar) <- list(name,c("lo","mean","up"))

xbar<- xbar*-1  #   get liberals on left -- estimates flipped


## party affiliation
name <- scan(file="C:/Documents and Settings/clinton/My Documents/Published Papers/CJR_PS2004/name.asc",
                   what=list(name=""),
                   sep="\n")$name
n <- length(name)
getparty <- function(x){
  indx <- regexpr("-",x)
  z <- substring(x,indx-1,indx-1)
  z
}
party <- getparty(name)
col <- rep(NA,n)
col[party=="D"] <- "blue"
col[party=="R"] <- "red"
col[party=="I"] <- "black"


###############################################################################################
#   Point Estimates and 95 % CI for posterior Means
#   Just National Journal Votes, No Bush
#
#   Similar to Figure 1 in Text

pdf(file="C:/Figure1.pdf",
    height=8,
    width=6)
par(mfrow=c(1,2),mar=c(2,4.1,2,.05))
indx <- order(xbar[,2])
plot(x=c(-6,0),y=c(1,50),type="n",axes=F,xlab="",ylab="",ylim=c(.5,50.5),yaxs="i")
abline(v=seq(-3,0,by=.5),lty=2,col=1,lwd=.5)
for(i in 1:50){
  lines(y=c(i,i),x=c(xbar[indx[i],"lo"],xbar[indx[i],"up"]),lwd=2.5,col=col[indx[i]])
  points(y=i,x=xbar[indx[i],"mean"],
         pch=183,cex=3,
         col=col[indx[i]])
  text(-4,i,dimnames(xbar)[[1]][indx[i]],cex=.5,adj=1)
}
axis(1,
     at=c(seq(-3,.5,by=1)),
     labels=c("-3","-2","-1","0"),
     cex=.85,lwd=2)
axis(3,
     at=c(seq(-3,.5,by=1)),
     labels=c("-3","-2","-1","0"),
     cex=.85,lwd=2)

plot(x=c(0,6),y=c(1,50),type="n",axes=F,xlab="",ylab="",ylim=c(.5,50.5),yaxs="i")
abline(v=seq(0,3,by=.5),lty=2,col=1,lwd=.5)
for(j in 1:50){
  i <- 50 + j
  lines(y=c(j,j),x=c(xbar[indx[i],"lo"],xbar[indx[i],"up"]),lwd=2.5,col=col[indx[i]])
  points(y=j,x=xbar[indx[i],"mean"],
         pch=183,cex=3,
         col=col[indx[i]])
  text(5,j,dimnames(xbar)[[1]][indx[i]],cex=.5,adj=1)
}
axis(1,
     at=c(seq(0,3,by=1)),
     labels=c("0","1","2","3"),
     cex=.85,lwd=2)
axis(3,
     at=c(seq(0,3,by=1)),
     labels=c("0","1","2","3"),
     cex=.85,lwd=2)
dev.off()


###############################################################################################
#   Get rank ordering
###############################################################################################

rankp<-matrix(NA,nrow=dim(NJ.key)[1],ncol=dim(NJ.key)[2])
for(i in 1:dim(NJ.key)[1]){
    rankp[i,]<-rank(-1*NJ.key[i,])
}


rankbar <- apply(rankp,2,function(x)quantile(x,c(.025,.50,.975)))
rankbar <- t(rankbar)
rankbar[,2] <- apply(rankp,2,mean)


dimnames(rankbar) <- list(name,c("lo","mean","up"))

#   Numbers are slightly different due to difference MCMC estimates but substantively identical

table(rankp[,60]==1)   # Prob Kerry most liberal
table(rankp[,81]==1)   # Prob Reed most liberal
table(rankp[,86]==1)   # Prob Sarbanes most liberal

table(rankp[,39 ] > rankp[,60])   # Prob Kerry more liberal than Edwards
table(rankp[,67 ] > rankp[,60])   # Prob Kerry more liberal than Lieberman
table(rankp[,67 ] > rankp[,39])   # Prob Edwards more liberal than Lieberman

table(rankp[,60 ] > rankp[,81])   # Prob Kerry more cons. than Reed
table(rankp[,60 ] > rankp[,86])   # Prob Kerry more cons than Sarbanes


###############################################################################################
#   Point Estimates and 95 % CI for Rank Ordering
#   Just National Journal Votes, No Bush
#
#   Similar to Figure 2 in Text

pdf(file="C:/Figure2.pdf",
    height=8,
    width=6)    
par(mfrow=c(1,2),mar=c(2,4.1,2,.05))
indx <- order(rankbar[,2])
plot(x=c(-25,55),y=c(1,50),type="n",axes=F,xlab="",ylab="",ylim=c(.5,50.5),yaxs="i")
abline(v=seq(1,55,by=10),lty=2,col=1,lwd=.5)
for(i in 1:50){
  lines(y=c(i,i),x=c(rankbar[indx[i],"lo"],rankbar[indx[i],"up"]),lwd=2.5,col=col[indx[i]])
  points(y=i,x=rankbar[indx[i],"mean"],
         pch=183,cex=3,
         col=col[indx[i]])
  text(0,i,dimnames(rankbar)[[1]][indx[i]],cex=.5,adj=1)
}
axis(1,
     at=c(seq(0,50,by=10)),
     labels=c("0","10","20","30","40","50"),
     cex=.85,lwd=2)
axis(3,
     at=c(seq(0,50,by=10)),
     labels=c("0","10","20","30","40","50"),
     cex=.85,lwd=2)

plot(x=c(50,140),y=c(1,50),type="n",axes=F,xlab="",ylab="",ylim=c(.5,50.5),yaxs="i")
abline(v=seq(50,100,by=10),lty=2,col=1,lwd=.5)
for(j in 1:50){
  i <- 50 + j
  lines(y=c(j,j),x=c(rankbar[indx[i],"lo"],rankbar[indx[i],"up"]),lwd=2.5,col=col[indx[i]])
  points(y=j,x=rankbar[indx[i],"mean"],
         pch=183,cex=3,
         col=col[indx[i]])
  text(123,j,dimnames(rankbar)[[1]][indx[i]],cex=.5,adj=1)
}
axis(1,
     at=c(seq(50,100,by=10)),
     labels=c("50","60","70","80","90","100"),
     cex=.85,lwd=2)
axis(3,
     at=c(seq(50,100,by=10)),
     labels=c("50","60","70","80","90","100"),
     cex=.85,lwd=2)
dev.off()

###############################################################################################
###############################################################################################
#  National Journal Votes Only
#  Include  Bush
###############################################################################################
###############################################################################################

idealout <- scan(file="C:/Documents and Settings/clinton/My Documents/ideal/Nat Journal/KeyBush/src/x.out")

NJ.key <- matrix(idealout,ncol=102,byrow=T)
NJ.key <- NJ.key[-seq(1:100),-1]      #   drop start values + burn-in (200,000), iter count

xbar <- apply(NJ.key,2,function(x)quantile(x,c(.025,.50,.975)))
xbar <- t(xbar)
xbar[,2] <- apply(NJ.key,2,mean)



name <- scan(file="C:/Documents and Settings/clinton/My Documents/Published Papers/CJR_PS2004/name.asc",
                   what=list(name=""),
                   sep="\n")$name

name<-c("BUSH",name)
                   
dimnames(xbar) <- list(name,c("lo","mean","up"))

#xbar<- xbar*-1  #   get liberals on left (if needed)


## party affiliation
name <- scan(file="C:/Documents and Settings/clinton/My Documents/Published Papers/CJR_PS2004/name.asc",
                   what=list(name=""),
                   sep="\n")$name
n <- length(name)
getparty <- function(x){
  indx <- regexpr("-",x)
  z <- substring(x,indx-1,indx-1)
  z
}
party <- getparty(name)

party<-c("R",party)

col <- rep(NA,n)
col[party=="D"] <- "blue"
col[party=="R"] <- "red"
col[party=="I"] <- "black"



###############################################################################################
#   Point Estimates and 95 % CI for posterior Means
#   Just National Journal Votes, include Bush
#
#   Basis for Figure 3 in Text

pdf(file="C:/Figure3basis.pdf",
    height=8,
    width=6) 
par(mfrow=c(1,2),mar=c(2,4.1,2,.05))
indx <- order(xbar[,2])
plot(x=c(-5,.5),y=c(1,51),type="n",axes=F,xlab="",ylab="",ylim=c(.5,51.5),yaxs="i")
abline(v=seq(-3,.5,by=.5),lty=2,col=1,lwd=.5)
for(i in 1:50){
  lines(y=c(i,i),x=c(xbar[indx[i],"lo"],xbar[indx[i],"up"]),lwd=2.5,col=col[indx[i]])
  points(y=i,x=xbar[indx[i],"mean"],
         pch=183,cex=3,
         col=col[indx[i]])
  text(-3.5,i,dimnames(xbar)[[1]][indx[i]],cex=.5,adj=1)
}
axis(1,
     at=c(seq(-3,.5,by=1)),
     labels=c("-3","-2","-1","0"),
     cex=.85,lwd=2)
axis(3,
     at=c(seq(-3,.5,by=1)),
     labels=c("-3","-2","-1","0"),
     cex=.85,lwd=2)

plot(x=c(-.5,5),y=c(1,51),type="n",axes=F,xlab="",ylab="",ylim=c(.5,51.5),yaxs="i")
abline(v=seq(-.5,3,by=.5),lty=2,col=1,lwd=.5)
for(j in 1:51){
  i <- 50 + j
  lines(y=c(j,j),x=c(xbar[indx[i],"lo"],xbar[indx[i],"up"]),lwd=2.5,col=col[indx[i]])
  points(y=j,x=xbar[indx[i],"mean"],
         pch=183,cex=3,
         col=col[indx[i]])
  text(4,j,dimnames(xbar)[[1]][indx[i]],cex=.5,adj=1)
}
axis(1,
     at=c(seq(0,3,by=1)),
     labels=c("0","1","2","3"),
     cex=.85,lwd=2)
axis(3,
     at=c(seq(0,3,by=1)),
     labels=c("0","1","2","3"),
     cex=.85,lwd=2)
dev.off()



###############################################################################################
#   Get rank ordering
###############################################################################################

rankp<-matrix(NA,nrow=dim(NJ.key)[1],ncol=dim(NJ.key)[2])
for(i in 1:dim(NJ.key)[1]){
    rankp[i,]<-rank(NJ.key[i,])
}

table(rankp[,1]==101)   # Prob Bush most conservative
table(rankp[,61]==1)   # Prob Kerry most liberal


rankbar <- apply(rankp,2,function(x)quantile(x,c(.025,.50,.975)))
rankbar <- t(rankbar)
rankbar[,2] <- apply(rankp,2,mean)

name <- scan(file="C:/Documents and Settings/clinton/My Documents/Published Papers/CJR_PS2004/name.asc",
                   what=list(name=""),
                   sep="\n")$name
name<-c("BUSH",name)

dimnames(rankbar) <- list(name,c("lo","mean","up"))


###############################################################################################
#   Point Estimates and 95 % CI for Rank Ordering
#   Just National Journal Votes, include Bush
#
#   Similar to Figure 3 in Text

pdf(file="C:/Figure3.pdf",
    height=8,
    width=6)
par(mfrow=c(1,2),mar=c(2,4.1,2,.05))
indx <- order(rankbar[,2])
plot(x=c(-25,55),y=c(1,51),type="n",axes=F,xlab="",ylab="",ylim=c(.5,51.5),yaxs="i")
abline(v=seq(1,55,by=10),lty=2,col=1,lwd=.5)
for(i in 1:50){
  lines(y=c(i,i),x=c(rankbar[indx[i],"lo"],rankbar[indx[i],"up"]),lwd=2.5,col=col[indx[i]])
  points(y=i,x=rankbar[indx[i],"mean"],
         pch=183,cex=3,
         col=col[indx[i]])
  text(-3,i,dimnames(rankbar)[[1]][indx[i]],cex=.5,adj=1)
}
axis(1,
     at=c(seq(0,50,by=10)),
     labels=c("0","10","20","30","40","50"),
     cex=.85,lwd=2)
axis(3,
     at=c(seq(0,50,by=10)),
     labels=c("0","10","20","30","40","50"),
     cex=.85,lwd=2)

plot(x=c(50,140),y=c(1,51),type="n",axes=F,xlab="",ylab="",ylim=c(.5,51.5),yaxs="i")
abline(v=seq(50,100,by=10),lty=2,col=1,lwd=.5)
for(j in 1:51){
  i <- 50 + j
  lines(y=c(j,j),x=c(rankbar[indx[i],"lo"],rankbar[indx[i],"up"]),lwd=2.5,col=col[indx[i]])
  points(y=j,x=rankbar[indx[i],"mean"],
         pch=183,cex=3,
         col=col[indx[i]])
  text(125,j,dimnames(rankbar)[[1]][indx[i]],cex=.5,adj=1)
}
axis(1,
     at=c(seq(50,100,by=10)),
     labels=c("50","60","70","80","90","100"),
     cex=.85,lwd=2)
axis(3,
     at=c(seq(50,100,by=10)),
     labels=c("50","60","70","80","90","100"),
     cex=.85,lwd=2)
dev.off()


###############################################################################################
###############################################################################################
#  All Sen 107 Votes 
#  Include Bush
###############################################################################################
###############################################################################################

readfunc <- function(file="data",first,last){
  foo <- scan(file=file,what=list(x=""),sep="\n")$x
  z <- substring(foo,first=first,last=last)
  z
}

readvotes <- function(file="data",first){
  foo <- readLines(file)
  n <- length(foo)
  m <- length(first:nchar(foo[1]))
  votes <- matrix(NA,n,m)
  for (i in 1:n){
    votes[i,] <- as.numeric(unlist(strsplit(substring(foo[i],first=first),"")))
  }
  votes
}

state<-as.numeric(readfunc(file="C:/Documents and Settings/clinton/My Documents/Published Papers/CJR_PS2004/sen107kh.ord.txt",9,11))
party<-as.numeric(readfunc(file="C:/Documents and Settings/clinton/My Documents/Published Papers/CJR_PS2004/sen107kh.ord.txt",21,23))
names<-readfunc(file="C:/Documents and Settings/clinton/My Documents/Published Papers/CJR_PS2004/sen107kh.ord",26,36)

idealout <- scan(file="C:/Documents and Settings/clinton/My Documents/ideal/Nat Journal/s107/src/x.out")

all.key <- matrix(idealout,ncol=103,byrow=TRUE)
all.key <- all.key[-seq(1:200),-1]      #   drop start values, iter count

xbar <- apply(all.key,2,function(x)quantile(x,c(.025,.50,.975)))
xbar <- t(xbar)
xbar[,2] <- apply(all.key,2,mean)

names<-names[-48]   #   drop barkley
party<-party[-48]   #   drop barkley

#   get names from Combine.Data.Rl
                   
dimnames(xbar) <- list(names,c("lo","mean","up"))

col <- rep(NA,dim(xbar)[1])
col[party==100] <- "blue"
col[party==200] <- "red"
col[party==328] <- "black"


###############################################################################################
#   Point Estimates and 95 % CI for posterior Means
#   Just National Journal Votes


par(mfrow=c(1,2),mar=c(2,4.1,2,.05))
indx <- order(xbar[,2])
plot(x=c(-1,0),y=c(1,51),type="n",axes=F,xlab="",ylab="",ylim=c(.5,51.5),yaxs="i")
abline(v=seq(-.5,0,by=.25),lty=2,col=1,lwd=.5)
for(i in 1:51){
  lines(y=c(i,i),x=c(xbar[indx[i],"lo"],xbar[indx[i],"up"]),lwd=2.5,col=col[indx[i]])
  points(y=i,x=xbar[indx[i],"mean"],
         pch=183,cex=3,
         col=col[indx[i]])
  text(-.5,i,dimnames(xbar)[[1]][indx[i]],cex=.65,adj=1)
}
axis(1,
     at=c(seq(-.5,.0,by=.25)),
     labels=c("-.5","-.25","0"),
     cex=.85,lwd=2)
axis(3,
     at=c(seq(-.5,.0,by=.25)),
     labels=c("-.5","-.25","0"),
     cex=.85,lwd=2)

plot(x=c(-.5,.75),y=c(1,51),type="n",axes=F,xlab="",ylab="",ylim=c(.5,51.5),yaxs="i")
abline(v=seq(0,.75,by=.25),lty=2,col=1,lwd=.5)
for(j in 1:51){
  i <- 51 + j
  lines(y=c(j,j),x=c(xbar[indx[i],"lo"],xbar[indx[i],"up"]),lwd=2.5,col=col[indx[i]])
  points(y=j,x=xbar[indx[i],"mean"],
         pch=183,cex=3,
         col=col[indx[i]])
  text(0,j,dimnames(xbar)[[1]][indx[i]],cex=.65,adj=1)
}
axis(1,
     at=c(seq(0,.5,by=.25)),
     labels=c("0",".25",".5"),
     cex=.85,lwd=2)
axis(3,
     at=c(seq(0,.5,by=.25)),
     labels=c("0",".25",".5"),
     cex=.85,lwd=2)
dev.off()

###############################################################################################
#   Get rank ordering
###############################################################################################

rankp<-matrix(NA,nrow=dim(all.key)[1],ncol=dim(all.key)[2])
for(i in 1:dim(all.key)[1]){
    rankp[i,]<-rank(all.key[i,])
}

rankbar <- apply(rankp,2,function(x)quantile(x,c(.025,.50,.975)))
rankbar <- t(rankbar)
rankbar[,2] <- apply(rankp,2,mean)

dimnames(rankbar) <- list(names,c("lo","mean","up"))

table(rankp[,1]==102)   # Prob Bush most conservative
table(rankp[,43]==1)   # Prob Kerry most liberal
rankbar[1,]
rankbar[43,]


table(rankp[,81]==1)   # Prob Reed most liberal
table(rankp[,86]==1)   # Prob Sarbanes most liberal

table(rankp[,39 ] > rankp[,60])   # Prob Kerry more liberal than Edwards
table(rankp[,67 ] > rankp[,60])   # Prob Kerry more liberal than Lieberman
table(rankp[,67 ] > rankp[,39])   # Prob Edwards more liberal than Lieberman

table(rankp[,60 ] > rankp[,81])   # Prob Kerry more cons. than Reed
table(rankp[,60 ] > rankp[,86])   # Prob Kerry more cons than Sarbanes



###############################################################################################
#   Point Estimates and 95 % CI for Rank Ordering
#   Just National Journal Votes, No Bush
#
#   Figure 4

pdf(file="C:/Figure4.pdf",
    height=8,
    width=6)
    
par(mfrow=c(1,2),mar=c(2,4.1,2,.05))
indx <- order(rankbar[,2])
plot(x=c(-25,55),y=c(1,51),type="n",axes=F,xlab="",ylab="",ylim=c(.5,51.5),yaxs="i")
abline(v=seq(1,55,by=10),lty=2,col=1,lwd=.5)
for(i in 1:51){
  lines(y=c(i,i),x=c(rankbar[indx[i],"lo"],rankbar[indx[i],"up"]),lwd=2.5,col=col[indx[i]])
  points(y=i,x=rankbar[indx[i],"mean"],
         pch=183,cex=3,
         col=col[indx[i]])
  text(-3,i,dimnames(rankbar)[[1]][indx[i]],cex=.5,adj=1)
}
axis(1,
     at=c(seq(0,50,by=10)),
     labels=c("0","10","20","30","40","50"),
     cex=.85,lwd=2)
axis(3,
     at=c(seq(0,50,by=10)),
     labels=c("0","10","20","30","40","50"),
     cex=.85,lwd=2)

plot(x=c(50,130),y=c(1,51),type="n",axes=F,xlab="",ylab="",ylim=c(.5,51.5),yaxs="i")
abline(v=seq(50,100,by=10),lty=2,col=1,lwd=.5)
for(j in 1:51){
  i <- 51 + j
  lines(y=c(j,j),x=c(rankbar[indx[i],"lo"],rankbar[indx[i],"up"]),lwd=2.5,col=col[indx[i]])
  points(y=j,x=rankbar[indx[i],"mean"],
         pch=183,cex=3,
         col=col[indx[i]])
  text(122,j,dimnames(rankbar)[[1]][indx[i]],cex=.5,adj=1)
}
axis(1,
     at=c(seq(50,100,by=10)),
     labels=c("50","60","70","80","90","100"),
     cex=.85,lwd=2)
axis(3,
     at=c(seq(50,100,by=10)),
     labels=c("50","60","70","80","90","100"),
     cex=.85,lwd=2)
dev.off()