pooling_termometro.R

## Please note, this script is entirely derivative of work done
## by Simon Jackman for his article `Pooling the Polls...'
## Aust. J. Pol. Science, 40(4): 499-517

## Libraries 
library(reshape)
library(rjags)
library(R2WinBUGS)
library(ggplot2)
library(scales)

#################################################################
## JAGS MODEL!
#################################################################

model1 <- function() {
 ## measurement model
 for(i in 1:nPolls){
       mu[i] <- house[org[i]] + alpha[date[i]]
       y[i] ~ dnorm(mu[i],prec[i])%_%T(0,1)
 }

 ## transition model (aka random walk prior)
 for(i in 2:nPeriods){
        mu.alpha[i] <- alpha[i-1]
        alpha[i] ~ dnorm(mu.alpha[i],tau)
 }

 ## priors
 tau <- 1/pow(sigma,2)   ## deterministic transform to precision
 sigma ~ dunif(0,.01)    ## uniform prior on standard deviation

 alpha[1] ~ dunif(0,.5) ## initialization of daily track

 for(i in 1:nOrgs){          ## vague normal priors for house effects
  house.star[i] ~ dnorm(0,.01)
	house[i] <- house.star[i] - mean(house.star)
 }

}
write.model(model1,con="model1.bug")

#################################################################
#################################################################

## Deal with locales
oldloc <- Sys.getlocale("LC_ALL")
if (!grepl("it_IT.utf8",oldloc)) {
	new <- Sys.setlocale(category="LC_ALL",locale="it_IT.utf8") 
}

## Read in data
dat <- read.delim2(file="termometro_polling.csv",sep=";",header=T,na.string=c("","NA","Null"),dec=",")
names(dat) <- c("Party","Day","Month","Company","Year","y")

##
dat$Party[ grepl("FLI",dat$Party) ] <- "FLI"

## Create date object
dat$Julian <- julian(as.Date(paste(dat$Day,tolower(dat$Month),dat$Year,sep="-"),
	"%d-%B-%Y"),
	origin=as.Date("2012-07-01"))

## Subset and order
dat <- dat[which(dat$Julian > 0),]
dat <- dat[!is.na(dat$Julian),]
dat <- dat[order(dat$Julian),]

## Reshape
dat <- dat[,c("Julian","Company","Party","y")]
dat <- cast(dat,Julian+Company~Party)
dat$Company <- factor(as.character(dat$Company))
the.parties <- c("PDL","LEGA",
	"UDC","FLI",
	"PD","SEL","IDV","FEDSIN",
	"M5S","SCMONTI")
nRespondents <- 1000 ## assumed

for (the.party in the.parties) {
	y <- dat[,the.party]/100
	var <- y*(1-y)/(nRespondents)

	forJags <- list(y=y[!is.na(y)],
		        prec=1/var[!is.na(y)],
		        date=(dat$Julian - min(dat$Julian) + 1),
		        org=as.integer(dat$Company),
				nOrgs = length(unique(dat$Company)),
		        nPolls=length(na.omit(y)),
		        nPeriods=length(min(dat$Julian):max(dat$Julian)),
		        alpha=c(rep(NA,
		          length(min(dat$Julian):max(dat$Julian))))
		        )

	initfunc <- function(){
	  house.star <- rnorm(n=forJags$nOrgs,sd=.075)
		nPeriods = length(min(dat$Julian):max(dat$Julian))
	 # alpha <- c(runif(n=1,min(y),max(y)),
	#	         runif(n=(nPeriods-1)))
	  sigma <- runif(n=1,0,.01)
	  list(house.star=house.star,
	#	   alpha=alpha,
		   sigma=sigma)
	}

	my.model <- jags.model("model1.bug",
		data=forJags,inits=initfunc,n.chains=1,
		n.adapt=100)
	out <- coda.samples(my.model,n.burnin=.2e6,
		variable.names=c("alpha","house","sigma"),
		n.iter=1e6,thin=1000)
	#	n.iter=1.2e6,
	#	n.thin=1000)
	#out <- jags.parfit(cl,
	#	data=forJags,
	#	params=c("alpha","house","sigma"),
	#	model=model1,
	#	inits = initfunc,	
	#	n.chains=2,
	#	n.adapt=100,
	#	n.burnin=.2e6,
	#	n.iter=1.2e6,
	#	n.thin=1000)

	filename = paste("jags_output_",the.party,".RData",sep="")
	save.image(file=filename)
	holder <- summary(out)

	## Create house effects -- want value, Hi, Lo
	house.fx.bar <- holder$statistics[grep("house",rownames(holder$statistics)),1]
	house.fx.lo <- holder$quantiles[grep("house",rownames(holder$quantiles)),1]
	house.fx.hi <- holder$quantiles[grep("house",rownames(holder$quantiles)),5]
	house.fx.df <- data.frame(Party=the.party,Organisation=levels(dat$Company),
		Lo=house.fx.lo,Mean=house.fx.bar,Hi=house.fx.hi)

	filename = paste("house_fx_",the.party,".csv",sep="")
	write.csv(house.fx.df,filename)

	my.alpha <- holder$statistics[grep("alpha",rownames(holder$quantiles)),1]
	my.dates <- seq(min(dat$Julian):max(dat$Julian))

	## Optionally plot
	plot(my.dates,my.alpha,type="p",pch=".",cex=2)
	points(forJags$date,y,type="p",pch=21,cex=1,col="blue")

	## Optionally plot house FX
	house.fx.df <- house.fx.df[order(house.fx.df$Mean),]
	house.fx.df$Organisation <- factor(house.fx.df$Organisation,
		levels=house.fx.df$Organisation,
		ordered=TRUE)

	filename <- paste("house_effects_",the.party,".pdf",sep="")
	pdf(file=filename)
	print(ggplot(house.fx.df,aes(x=Organisation,y=Mean,ymin=Lo,ymax=Hi)) + 
		scale_y_continuous("Mean house effect",labels=percent) + 
		geom_pointrange() + geom_hline(yintercept=0) + theme_bw() + opts(title=the.party,axis.text.x=theme_text(angle=90)))
	dev.off()
}