library(ggplot2)
library(plyr)
library(Matrix)
library(BradleyTerry2)

## Read in the data
ranks<-read.delim("rankings.csv",header=T,as.is=TRUE,sep="\t")
ranks$X<-NULL

## Add ranks for years in which we have just scores
## 
ranks<-ddply(ranks,.(Year),function(df) {
	if (all(is.na(df$Lucid))) {
		df$Lucid<-rank(df$LucidScore)
		df$Star<-rank(df$StarScore)
		df$Domestic<-rank(df$DomesticScore)
		df$PolRank<-rank(df$PolScore)
		df
	}
	else {
		df
	}
})

## Show the year by year correlation of 
## Economic and political ranking
pdf(file="pol-econ-cor.pdf")
par(mfrow=c(2,3))
## For each year
for (i in unique(ranks$Year)) {
	if (any(!is.na(ranks$EconRank[ranks$Year==i]))) {
		plot(ranks$EconRank[ranks$Year==i],ranks$PolRank[ranks$Year==i],
			pch=19,ylim=c(0,20),xlim=c(0,20),
			ylab="Political ranking",xlab="Economic ranking",
			main=paste("Year: ",i))
		abline(lm(ranks$EconRank[ranks$Year==i]~ranks$PolRank[ranks$Year==i]),lty=2,col="red")
		my.cor<-cor(ranks$PolRank[ranks$Year==i],ranks$EconRank[ranks$Year==i],method="spearman",use="pairwise")
		text(19,19,paste("r=",round(my.cor,digits=3)))
	}
}
## And now for all years
plot(jitter(ranks$EconRank),jitter(ranks$PolRank),pch=19,ylim=c(0,20),xlim=c(0,20),ylab="Political ranking",xlab="Economic ranking",main="All years")
abline(lm(ranks$EconRank~ranks$PolRank),lty=2,col="red")
my.cor<-cor(ranks$PolRank,ranks$EconRank,method="spearman",use="pairwise")
text(19,19,paste("r=",round(my.cor,digits=3)))

dev.off()


## Now generate pairwise contests
pairwise.df<-expand.grid(Player1=unique(ranks$Minister),Player2=unique(ranks$Minister),Year=unique(ranks$Year))
pairwise.df$Wins1<-pairwise.df$Wins2<-pairwise.df$Tied<-NA
## Eliminate contests-against-self
pairwise.df<-pairwise.df[which(pairwise.df$Player1 != pairwise.df$Player2),]


for (i in 1:nrow(pairwise.df)) {

	p1.ranks<-ranks$PolRank[which(ranks$Minister==pairwise.df$Player1[i] & ranks$Year == pairwise.df$Year[i])]
	p2.ranks<-ranks$PolRank[which(ranks$Minister==pairwise.df$Player2[i] & ranks$Year == pairwise.df$Year[i])]
	## Remember, higher ranks are bad, so the usual pairwise comparison is arse-backwards
	if (length(p1.ranks) & length(p2.ranks)) {
		p1btp2<-p1.ranks < p2.ranks
		p2btp1<-p1.ranks > p2.ranks
		tied <- p1.ranks == p2.ranks
		pairwise.df$Wins1[i]<- sum(p1btp2,na.rm=T)
		pairwise.df$Wins2[i]<- sum(p2btp1,na.rm=T)
		pairwise.df$Tied[i] <- sum(tied)
	} else {
		pairwise.df$Wins1[i]<- NA
		pairwise.df$Wins2[i]<- NA
		pairwise.df$Tied[i]<- NA
	}
}

## Remove extraneous rows
pairwise.df<-pairwise.df[!is.na(pairwise.df$Wins1),]
pairwise.df<-pairwise.df[!is.na(pairwise.df$Wins2),]

## Account for ties?
pairwise.df$Wins1.adj<-pairwise.df$Wins1 + pairwise.df$Tied/2
pairwise.df$Wins2.adj<-pairwise.df$Wins2 + pairwise.df$Tied/2


contest<-list(winner=pairwise.df)

baseline.BTm <- BTm(cbind(Wins1, Wins2), 
	Player1, Player2, 
	formula=~..,
	data = contest,family=binomial(probit),refcat="Peer Steinbrück",br=TRUE)



## Now do the Bradley-Terry model
predictors<-read.delim("invariant_minister_data.csv",header=T,sep="\t")

plot.abilities<-function(x) { 
	BTm.summary<-as.data.frame(BTabilities(x))
	BTm.summary$Hi<-BTm.summary$"ability" + 1.645 * BTm.summary$"s.e."
	BTm.summary$Lo<-BTm.summary$"ability" - 1.645 * BTm.summary$"s.e."
	BTm.summary$Minister<-rownames(BTm.summary)
	BTm.summary$Minister<-gsub("\\.\\.","",as.character(BTm.summary$Minister))
	BTm.summary$Minister<-gsub("Player","",as.character(BTm.summary$Minister))
	BTm.summary<-BTm.summary[order(BTm.summary$"s.e.",decreasing=TRUE),]

	BTm.summary$Minister<-factor(BTm.summary$Minister,levels=BTm.summary$Minister,ordered=TRUE)
	BTm.summary$Minister<-reorder(BTm.summary$Minister,BTm.summary$ability)
	BTm.summary<-merge(BTm.summary,predictors,by="Minister",all.x=T,all.y=F)
	my.plot <- ggplot(subset(BTm.summary,BTm.summary$"s.e."<100),aes(x=Minister,y=ability,ymax=Hi,ymin=Lo,color=iso3c)) + geom_pointrange() + coord_flip() + theme_bw()
	return(my.plot)
}

baseline.plot<-plot.abilities(baseline.BTm)

pdf(file="baseline_abilities.pdf")
print(baseline.plot)
dev.off()
	BTm.summary<-as.data.frame(BTabilities(baseline.BTm))
	BTm.summary$Hi<-BTm.summary$"ability" + 1.645 * BTm.summary$"s.e."
	BTm.summary$Lo<-BTm.summary$"ability" - 1.645 * BTm.summary$"s.e."
	BTm.summary$Minister<-rownames(BTm.summary)
	BTm.summary$Minister<-gsub("\\.\\.","",as.character(BTm.summary$Minister))
	BTm.summary$Minister<-gsub("Player","",as.character(BTm.summary$Minister))
	BTm.summary<-BTm.summary[order(BTm.summary$"s.e.",decreasing=TRUE),]

	BTm.summary$Minister<-factor(BTm.summary$Minister,levels=BTm.summary$Minister,ordered=TRUE)
	BTm.summary$Minister<-reorder(BTm.summary$Minister,BTm.summary$ability)
	BTm.summary<-merge(BTm.summary,predictors,by="Minister",all.x=T,all.y=F)
BTm.summary$iso3c<-factor(BTm.summary$iso3c,levels=unique(BTm.summary$iso3c),ordered=TRUE)
BTm.summary$iso3c<-reorder(BTm.summary$iso3c,BTm.summary$ability)
country.plot <- ggplot(subset(BTm.summary,BTm.summary$"s.e."<100),aes(x=iso3c,y=ability)) + geom_boxplot() + geom_point(position="jitter") + coord_flip() + theme_bw()

## Okay, format it the chameleons way

## Eliminate ties
winner.a<-as.character(pairwise.df[which(pairwise.df$Wins1==1),"Player1"])
loser.a<-as.character(pairwise.df[which(pairwise.df$Wins1==1),"Player2"])

winner.b<-as.character(pairwise.df[which(pairwise.df$Wins2==1),"Player2"])
loser.b<-as.character(pairwise.df[which(pairwise.df$Wins2==1),"Player1"])

winner<-data.frame(ID=c(winner.a,winner.b))
loser<-data.frame(ID=c(loser.a,loser.b))
my.levels<-unique(c(as.character(winner$ID),as.character(loser$ID)))
winner$ID<-factor(winner$ID,levels=my.levels)
loser$ID<-factor(loser$ID,levels=my.levels)

## Get predictors
predictors<-read.delim("invariant_minister_data.csv",header=T,sep="\t")
predictors$Minister<-factor(predictors$Minister,levels=levels(winner$ID))
## MISSING ENTRIES!
predictors<-predictors[!is.na(predictors$Minister),]
rownames(predictors)<-predictors$Minister

contest<-list(winner=winner,loser=loser,predictors=predictors)

summary(pred.BTm <- BTm(player1 = winner, player2 = loser,
       formula = ~  technocrat[ID] + ID , id = "ID",
	 data = contest))

foobar<-plot.abilities(pred.BTm)

summary(pred.BTm <- BTm(player1 = winner, player2 = loser,
       formula = ~  technocrat[ID] + female[ID] + (1|ID), id = "ID",family=binomial(link="probit"),
	 data = contest))