#########################################################
##  Josh Clinton
##  Princeton University
##
##  Example Code running IDEAL
##  Ohio State Brownbag: January 2008
#########################################################

#	Must install the ideal package from CRAN (with dependencies)
#	Must install the nominate package from CRAN (with dependencies)


require(pscl)				#	Load the pscl package
data(s109)								
is(s109,"rollcall")			# logical: is s109 a rollcall object?	
s109        				#   109th Senate   (built-in)
summary(s109)				#	Give summary stats of the 109th Senate
summary(s109,verbose=TRUE)	#	Give more extensive summary stats

sc9497      				#   1994-1997 Supreme Court   (built-in)


#   Read in Roll Call Files from Voteview website directly
#	(Need an active internet connect for this to work)
#	(Could also readKH from a file on your local computer by replacing the path)
h107 <- readKH("ftp://voteview.com/hou107kh.ord",desc="107th U.S. House of Representatives")

#	run: help(ideal) to get a list of options

id1 <- ideal(s109,			#	run MCMC Quadratic model on s109 data
    d=1,					#	1 dimensional
    normalize=TRUE,			#	Normalize the space s.t. x has mean 0, var. 1
    store.item=TRUE,		#	Store the alpha and betas for access later
    maxiter=500,			#	Run 500 iterations (way too short)
    burnin=100,				#	Discard the first 100 as "burnin" to the posterior mode
    thin=10,				#	Only analyze 1 out of every 10 posterior draws to explore posterior space
    verbose=TRUE)			#	Inform me what the code is doing at every step
    
summary(id1)				#	Summarize the output of the roll call analysis

names(id1)      #   See what objects are kept

id1$xbar        #   Posterior means (ideal points)
id1$x           #   Posterior draws
id1$betabar     #   Posterior means (item parameters) [beta,alpha]

cpmean<-id1$betabar[,2]/id1$betabar[,1]     #   Calculate cutpoint mean (Cauchy rv)
                                            #   alpha/beta

#   Characterize Congressional Voting Behavior

plot(density(id1$xbar[,1]),                 #   Plot distribution of x's and cp's
    main="109th Senate (2005-2006)",xlab="IDEAL Estimates")
rug(cpmean, ticksize=.1)					#	Add cutpoints along the bottom

#   Do a little nicer representation 

plot(density(id1$xbar[s109$legis$party=="D",1])	# Democrat distribution
    col="RED",lwd=2,xlim=c(-3,3),ylim=c(0,1.75),
    main="109th Senate (2005-2006)",xlab="IDEAL Estimates")

points(density(id1$xbar[s109$legis$party=="R",1])$x,  # Republican distribution
    density(id1$xbar[s109$legis$party=="R",1])$y,
    col="BLUE",lwd=2,type="l")

rug(cpmean)									        # plot cutpoint means
points(density(cpmean)$x,density(cpmean)$y,type="l") # plot distribution of cutpoint means


#   Informally "Assess" Convergence
#   Use idealtoMCMC for formal diagnostics

tracex(id1, legis="KENN", d=1, 
    conf.int=0.95,
    showAll = FALSE, 
    burnin=NULL,
    span=.25,
    legendLoc="topright")

tracex(id1, legis="FRIST", d=1, 
    conf.int=0.95,
    showAll = FALSE, 
    burnin=NULL,
    span=.25,
    legendLoc="topright")
    
tracex(id1,legis=c("KENN","BOX","KYL"))    

#   Summarize estimates for individual legislators

plot.ideal(id1, d=1,
    conf.int=0.95, 
    burnin=NULL,
    showAllNames = TRUE)

plot.ideal(id1, d=1,
    conf.int=0.95, 
    burnin=NULL,
    showAllNames = FALSE)

#	Now Get NOMINATE

library(wnominate)

help(wnominate)					#	look at the possible options

s109nom<-wnominate(s109, polarity=c(44,50))		# Kyl, Lieberman
plot(s109nom)

s109nom$legislators   # The objects relevant for legislators from nominate

s109nom$legislators$coord1D #	first dimension NOMINATE ideal points 

dev.off()
plot(id1$xbar,s109nom$legislators$coord1D)   # Plot NOINATE against IDEAL
cor(id1$xbar,s109nom$legislators$coord1D) 	 # Get correlation