### Week 9: Uncertainty
### Code from class

require(foreign)
require(Zelig)

## Examples from Homework 6
## titanic data qn3
titanic <- read.dta("titanic.dta")
levels(titanic$class) <- c("first","second","third","crew")
z.out <- zelig(survived ~ age+sex+class, model="logit", data=titanic)
summary(z.out)
x.kate <- setx(z.out, ageadults=1, sexman=0,
               classsecond=0, classthird=0, classcrew=0)
x.kate[1,] <- c(1,1,0,0,0,0)
x.leo <- setx(z.out, ageadults=1, sexman=1,
              classsecond=0, classthird=1, classcrew=0)
x.leo[1,] <- c(1,1,1,0,1,0)
summary(s.out <- sim(z.out, x=x.leo, x1=x.kate))

## economic_bills data qn4b
ecbills <- read.dta("economic_bills.dta")
z.out <- zelig(status ~ cabinet + vpdi_LH92economic + xland,
               model="logit", data=ecbills)
x.out <- setx(z.out)
x.out[1,] <- c(1,0,0,0,0,1)
summary(sim(z.out, x.out))
# for comparison:
log2 <- glm(status ~ cabinet + vpdi_LH92economic + xland, family=binomial,
            data=ecbills)
predict(log2,new=data.frame(cabinet=0,vpdi_LH92economic=0,xland="UK"),
        type="response", se=T)

## economic_bills data qn4c
x.out[1,] <- c(1,1,5,1,0,0)
summary(sim(z.out, x.out))
# for comparison:
predict(log2,new=data.frame(cabinet=1,vpdi_LH92economic=5,xland="FRA"),
        type="response",se=T)

## economic_bills data qn4d
(x.out <- setx(z.out, vpdi_LH92economic=seq(0,10.4,.1)))
x.out[,2] <- 0
x.out[,5] <- 1
s.out <- sim(z.out, x.out)
plot.ci(s.out)
lines(seq(0,10.4,.1), apply(s.out$qi$ev,2,mean))


## 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(1,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(1,1))


## 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")