###############################################################
# Replications for Figures 5 and 6 and Tables 1 and 2
# A Quantitative Method for Substantive Robustness Assessment
# by Justin Esarey and Nathan Danneman
#
# version: 5/11/2014
###############################################################


##################################
# Replicate Table 2
##################################


rm(list=ls())
library(foreign)
library(car)
library(snow)
setwd("C:/Dropbox/jesarey_documents/Substantive Robustness Project/Replication Materials")
source("ols_cstar_loss.r")
source("me_cstar_loss.r")
source("clusterse.r")

cl<-makeCluster(2, type="SOCK")
clusterExport(cl, "cstarme")
clusterExport(cl, "clx")

dat.in<-read.dta(file="legislative_cleaned.dta")
attach(dat.in)
dat<-na.omit(data.frame(enep1, proximity1, enpres,  proximity1_enpres, eneg, logmag, logmag_eneg, country))
dat$country<-factor(dat$country)
detach(dat.in)
attach(dat)
goldermod<-lm(enep1~proximity1+enpres+proximity1_enpres+eneg+logmag+logmag_eneg)
summary(goldermod)  ## This is Table 2
coef<-coefficients(goldermod)
vcv<-clx(goldermod, 1, country)
sqrt(diag(vcv))


##################################
# Replicate Figure 6
##################################


# plot derivative values
Z.fits<-seq(from=0.05, to=6, by=0.1)
# this is the main effect
dYdX<-coef[2]+Z.fits*coef[4]
plot(dYdX~Z.fits, type="l", ylim=c(-6, 4), xlab="# of Effective Presidential Candidates", ylab="ME of Presidential Elections", main=c("Marginal Effect of Proximate Presidential", "Elections on # of Legislative Parties"), lwd=2, col="black")
legend("bottomright", lty=c(1,2), legend=c("marginal effect", "95% CI"))

# what is the standard error of this quantity?
# exact calculation 
dYdX.se<-sqrt(vcv[2,2]+Z.fits^2*vcv[4,4]+2*(Z.fits)*vcv[2,4])
# 95% CI
lines(dYdX+1.96*dYdX.se~Z.fits, lty=2, lwd=2)
lines(dYdX-1.96*dYdX.se~Z.fits, lty=2, lwd=2)
# show the zero line
abline(h=0, lty=3)

dev.copy2eps(file="golderplot.eps")     ##  This is Figure 6 Panel (a)



# draw out of the asymptotic distribution of betas, fixed sigma (hat)
require(mvtnorm)
beta.draws<-rmvnorm(20000, mean=coef, sigma=vcv) # note: mean and vcv come out of regression results
dYdX.draws<-beta.draws[,2]+as.matrix(Z.fits)%*%beta.draws[,4]
dYdX.lo<-apply(dYdX.draws, MARGIN=1, FUN=quantile, probs=0.025)
dYdX.hi<-apply(dYdX.draws, MARGIN=1, FUN=quantile, probs=0.975)

r.test<-c(2, 4, 6, 10, 15, 20)
cstar.plot<-matrix(data=NA, nrow=length(r.test), ncol=60)
for(i in 1:length(r.test)){
  if(i==1){pb<-txtProgressBar(min=0, max=length(r.test), style=3)}
  setTxtProgressBar(pb, i)
  cstar.plot[i,]<-parApply(cl=cl, FUN=cstarme, MARGIN=1, X=dYdX.draws, r=r.test[i])
}
close(pb)


plot(dYdX~Z.fits, type="l", ylim=c(-4, 2), xlab="# of Effective Presidential Candidates", ylab="ME of Presidential Elections", main="Substantive Robustness Plot", lwd=2, col="black")

cstar.smooth<-matrix(data=NA, nrow=length(r.test), ncol=60)
for(i in 1:length(r.test)){
  lo.mod<-loess(cstar.plot[i,]~Z.fits, degree=1, span=0.4)
  cstar.smooth[i,]<-predict(lo.mod, newdata=Z.fits)
}
polygon(c(Z.fits, rev(Z.fits)), c(pmax(dYdX, 0), rev(pmin(dYdX,0))),
        col = paste("gray",50,sep=""), border = NA)

for(i in 1:6){
  lines(cstar.smooth[i,]~Z.fits, col=paste("gray",50+6*i,sep=""))
  polygon(c(Z.fits, rev(Z.fits)), c(pmax(cstar.smooth[i,], 0), rev(pmin(cstar.smooth[i,],0))),
          col = paste("gray",50+6*i,sep=""), border = NA)
  
}
lines(dYdX~Z.fits, col="black", lwd=2)

legtext<-list()
for(j in 1:6){legtext[[j]]<-substitute(gamma <= J, list(J=r.test[j]))}
legtext[[7]]<-substitute(gamma > J, list(J=r.test[6]))

legend("bottomright", legend=parse(text=legtext), fill=paste("gray",50+6*0:6,sep=""), text.width=0.6)

abline(h=0, lty=3)
text(0, 1.95, cex=0.8, pos=4, "* shading indicates ME is substantively significant")
text(0, 1.75, cex=0.8, pos=4, expression(paste("   for researcher with corresponding loss aversion ", gamma)))

dev.copy2eps(file="golder-robust.eps")    ## This is Figure 6 Panel (b)

stopCluster(cl)
detach(dat)






##################################
# Replicate Table 1
##################################

rm(list=ls())
library(foreign)
library(car)
library(snow)
source("ols_cstar_loss.r")

dat<-read.dta(file="clinton_repdata.dta")

dat.0<-subset(dat, subset=(pty==0))
dat.1<-subset(dat, subset=(pty==1))
mod.0<-lm(x106kmean~pctspid+pctnspid, weight=n, data=dat.0)
mod.1<-lm(x106kmean~pctspid+pctnspid, weight=n, data=dat.1)

##### This is Table 1
summary(mod.0)  # democratic legislators, key votes
summary(mod.1)  # republican leglslators, key votes



##################################
# Replicate Figure 5
##################################

r<-seq(from=1, to=6, by=0.1)
cstar.0<-matrix(data=NA, nrow=length(r), ncol=2)
cstar.1<-matrix(data=NA, nrow=length(r), ncol=2)
for(i in 1:length(r)){
  cstar.0[i,]<-cstar(mod.0, r=r[i])[2:3]
  cstar.1[i,]<-cstar(mod.1, r=r[i])[2:3]
}

par(mar=c(5,6,4,2))
plot(cstar.0[,1]~r, type="l", main=c("Estimated c* (effect of Democrat constituent", "ideology on Democratic legislators' key votes)"), xlab="level of loss aversion (gamma)", ylab="", ylim=c(-0.02, 2.22))
text(x=3.5, y=1.25, labels="act on null")
text(x=2, y=0.5, labels="act on alternative")
arrows(1.8, 0.4, 1.15, 0.05, length=0.15)
mtext(side=2, line=4, "minimum substantively meaningful relationship")
mtext(side=2, line=3, "d Democrat legislator ID / d Democrat constituent ID (c*)")


dev.copy2eps(file="clinton_DonD.eps")      ## This is Figure 5 Panel (d)


par(mar=c(5,6,4,2))
plot(cstar.0[,2]~r, type="l", main=c("Estimated c* (effect of Republican constituent", "ideology on Democrat legislators' key votes)"), xlab="level of loss aversion (gamma)", ylab="", ylim=c(0, 2.22))
text(x=4, y=2.1, labels="act on null")
text(x=3, y=1.6, labels="act on alternative")
mtext(side=2, line=4, "minimum substantively meaningful relationship")
mtext(side=2, line=3, "d Democrat legislator ID / d Republican constituent ID (c*)")


dev.copy2eps(file="clinton_RonD.eps")     ## This is Figure 5 Panel (c)



par(mar=c(5,6,4,2))
plot(cstar.1[,2]~r, type="l", main=c("Estimated c* (effect of Democrat constituent", "ideology on Republican legislators' key votes)"), xlab="level of loss aversion (gamma)", ylab="", ylim=c(0, 2.22))
text(x=4, y=0.6, labels="act on null")
text(x=2.5, y=0.14, labels="act on alternative")
mtext(side=2, line=4, "minimum substantively meaningful relationship")
mtext(side=2, line=3, "d Republican legislator ID / d Democrat constituent ID (c*)")


dev.copy2eps(file="clinton_DonR.eps")     ## This is Figure 5 Panel (b)


par(mar=c(5,6,4,2))
plot(cstar.1[,1]~r, type="l", main=c("Estimated c* (effect of Republican constituent", "ideology on Republican legislators' key votes)"), xlab="level of loss aversion (gamma)", ylab="", ylim=c(0, 2.22))
text(x=4, y=1.5, labels="act on null")
text(x=2.2, y=1.1, labels="act on alternative")
mtext(side=2, line=4, "minimum substantively meaningful relationship")
mtext(side=2, line=3, "d Republican legislator ID / d Republican constituent ID (c*)")


dev.copy2eps(file="clinton_RonR.eps")     ## This is Figure 5 Panel (a)



## This is where the substantive effects reported in the text
## of the Clinton replication are calculated

# effect of republican constituents on both parties
# first: democratic legislators
sd((subset(pctnspid, pty==0)))
Fn<-ecdf(subset(x106kmean, pty==0))
Fn(-.7687515+0.107887*1)-Fn(-.7687515)

# second: republican legislators
sd((subset(pctspid, pty==1)))
Fn<-ecdf(subset(x106kmean, pty==1))
Fn(.8055426 +0.07063387*1)-Fn(.8055426)


# effect of democrat constituents on both parties
# first: democrat legislators
sd((subset(pctspid, pty==0)))
Fn<-ecdf(subset(x106kmean, pty==0))
Fn(-.7687515+0.10315*0)-Fn(-.7687515)

# second: republican legislators
sd((subset(pctnspid, pty==1)))
Fn<-ecdf(subset(x106kmean, pty==1))
Fn(.8055426 +0.0859733*0.2)-Fn(.8055426)