#Script file based on code supplied by Ken Benoit and adjusted for our purposes.
#
library(effects)
    library(graphics)
    library(car)
    library(foreign)
   library(Zelig)


      dir<-setwd('d:\\classes\\pol290g\\Spring2009\\lectures')
      # dir<-setwd('c:\\braddocs\\classes\\pol290g\\Spring2009\\lectures') 

#Poisson Model using Zelig

## predicted values on Weede dataset from Benoit (1996)
weede <- read.dta("weede.dta")
z.out <- zelig(ssal6080 ~
               fh73+lpopln70+lmilwp70, model="poisson", data=weede)
(x.out <- setx(z.out, fh73=2:14))
s.out <- sim(z.out, x=x.out)
summary(s.out)

## replicate part of Table 3 from Benoit (1996)
z.tab2NBpoldem <- zelig(butterw ~ poldem65, model="negbin", data=weede)
x.tab2NBpoldem <- setx(z.tab2NBpoldem, poldem65=c(0,20,55,85,100))
s.tab2NBpoldem <- sim(z.tab2NBpoldem, x=x.tab2NBpoldem)
cbind(apply(s.tab2NBpoldem$qi$ev, 2, mean),
      apply(s.tab2NBpoldem$qi$ev, 2, sd))
z.tab2NBfh73 <- zelig(butterw ~ fh73, model="negbin", data=weede)
x.tab2NBfh73 <- setx(z.tab2NBfh73, fh73=c(2,4,7,12,14))
s.tab2NBfh73 <- sim(z.tab2NBfh73, x=x.tab2NBfh73)
cbind(apply(s.tab2NBfh73$qi$ev, 2, mean),
      apply(s.tab2NBfh73$qi$ev, 2, sd))

## replicate top part of Figure 1 from Benoit (1996)
x.tab2NBpoldem <- setx(z.tab2NBpoldem, poldem65=seq(0,100,1))
s.tab2NBpoldem <- sim(z.tab2NBpoldem, x=x.tab2NBpoldem)
x.tab2NBfh73 <- setx(z.tab2NBfh73, fh73=2:14)
s.tab2NBfh73 <- sim(z.tab2NBfh73, x=x.tab2NBfh73)
par(mfrow=c(2,2), mar=c(4,4,1,1))
plot.ci(s.tab2NBpoldem, xlab="POLDEM 1965", ylab="Butterworth Wars",
        ylim=c(0,7))
points(weede$poldem65, weede$butterw)
plot.ci(s.tab2NBfh73, xlab="Freedom House 1973", ylab="Butterworth Wars",
        ylim=c(0,7))
points(weede$fh73, weede$butterw)
par(mfrow=c(2,2))


## replicate Table 5 Benoit and Marsh (2009) PRQ
# note: convert.factors=F since this makes dummy vars 0/1 numeric
dail <- read.dta("dailprobit.dta", convert.factors=F)

z.out <- zelig(wonseat ~ pspend_total*incumb+m, model="probit", data=dail)

x.lo <- setx(z.out, pspend_total=5, incumb=0, m=4)
x.hi <- setx(z.out, pspend_total=15, incumb=0, m=4)
summary(s.out <- sim(z.out, x=x.lo, x1=x.hi))

x.lo <- setx(z.out, pspend_total=0, incumb=1, m=4)
x.hi <- setx(z.out, pspend_total=5, incumb=1, m=4)
summary(s.out <- sim(z.out, x=x.lo, x1=x.hi))

x.lo <- setx(z.out, pspend_total=5, incumb=1, m=4)
x.hi <- setx(z.out, pspend_total=10, incumb=1, m=4)
summary(s.out <- sim(z.out, x=x.lo, x1=x.hi))

x.lo <- setx(z.out, pspend_total=10, incumb=1, m=4)
x.hi <- setx(z.out, pspend_total=15, incumb=1, m=4)
summary(s.out <- sim(z.out, x=x.lo, x1=x.hi))

x.lo <- setx(z.out, pspend_total=5, incumb=1, m=4)
x.hi <- setx(z.out, pspend_total=15, incumb=1, m=4)
summary(s.out <- sim(z.out, x=x.lo, x1=x.hi))

## replicate Figure 2 Benoit and Marsh (2009) PRQ
dail <- read.dta("dailprobit.dta", convert.factors=F)
z.out <- zelig(wonseat ~ pspend_total*incumb+m, model="probit", data=dail)
x.incumb <- setx(z.out, pspend_total=seq(0,30,.5), incumb=1, m=4)
x.chall <- setx(z.out, pspend_total=seq(0,30,.5), incumb=0, m=4)
x.chall[1,5] <- .0001
s.out <- sim(z.out, x=x.incumb, x1=x.chall)
plot.ci(s.out, xlab="% Candidate Spending in Constituency",
        ylab="Probability of Winning a Seat")
text(5,.7,"Incumbents", col="red")
text(17,.4,"Challengers", col="blue")
abline(h=.5, lty="dashed", col="grey60")

## plot an ROC plot comparing challengers v. incumbent predictions
dail.incumb <- subset(dail, incumb==1, select=c(wonseat,pspend_total,incumb,m))
dail.chall <- subset(dail, incumb==0, select=c(wonseat,pspend_total,incumb,m))
z.out.i <- zelig(wonseat ~ pspend_total+m, model="probit", data=dail.incumb)
z.out.c <- zelig(wonseat ~ pspend_total+m, model="probit", data=dail.chall)
rocplot(z.out.i$y, z.out.c$y, fitted(z.out.i), fitted(z.out.c),
        lty1="solid", lty2="solid", col2="blue", col1="red")
text(.6,.55,"Incumbents",col="red")
text(.8,.85,"Challengers",col="blue")