############################################################### # Replications for Figures 1, 2, and 3 # A Quantitative Method for Substantive Robustness Assessment # by Justin Esarey and Nathan Danneman # # version: 5/11/2014 ############################################################### ################################## # Replicate Figure 1 ################################## rm(list=ls()) setwd("C:/Dropbox/jesarey_documents/Substantive Robustness Project/Replication Materials") set.seed(437823) x<-seq(from=-1, to=2, by=0.01) y<-dnorm(x,mean=0.25,sd=.25/.67) plot(y~x,type="l", xlab="Proportion Change in Cancer Probability After Exposure", ylab="Posterior Belief Density", main=c("Hypothetical Weight Loss Drug Analysis","Change in Cancer Probability After Treatment")) abline(v=0, lty=2) dev.copy2eps(file="cancerexample.eps") # This is Figure 1 Panel (a) x<-seq(from=-1, to=2, by=0.01) y<-dnorm(x,mean=0.25,sd=.25/.67) plot(y~x,type="l", xlab="Proportion Difference in Intrastate Conflict Probability", ylab="Posterior Belief Density", main=c("Hypothetical Conflict Analysis: Difference in","War Probability (Non-Democracy - Democracy)")) abline(v=0, lty=2) dev.copy2eps(file="conflictexample.eps") # This is Figure 1 Panel (b) ################################## # Replicate Figure 3 ################################## simplestati<-function(b,v,r){ # "immediate" form of simple kink-slope statistic # give it a beta and a standard error and an alpha, it returns a "loss point" t s<-sign(b) b<-abs(b) # reflect analysis into the positive domain (interpretation is the same) utility<-function(t,a,x,m,sigma){ k<-as.numeric((x-t)>0)*2-1 out<-a^(-k)*dnorm(x, mean=m, sd=sigma)*(x-t); return(out) } maximand<-function(t,a,m,sigma){as.numeric(integrate(utility, -Inf, Inf, t=t, a=a, m=m, sigma=sigma)[1])} out<-c() for(i in 1:length(b)){ if(floor(i/10)==(i/10)){print(c("Iteration ",i))} if(b[i]!=0){out[i]<- sign(s[i])*uniroot(maximand, interval=c(-100,100), a=r, m=b[i], sigma=v[i])$root} else{out[i]<- uniroot(maximand, interval=c(-100,100), a=r, m=b[i], sigma=v[i])$root} } return(out) } par(mar=c(5,6,4,2)) rvec<-seq(from=1, to=2.5, by=0.05) c_star<-pmax(100*unlist(lapply(rvec, FUN=simplestati,b=0.25, v=.25/.67)),0) plot(c_star~rvec, type="l", main=c("Estimated c* (effect of democracy on war probability)", "by level of loss aversion"), xlab="level of loss aversion (gamma)", ylab="") text(x=2, y=18, labels="act on null") text(x=1.4, y=5, labels="act on alternative") mtext(side=2, line=4, "minimum substantively meaningful relationship") mtext(side=2, line=3, "% change in war probability, non-democracy - democracy (c*)") dev.copy2eps(file="conflictrobust.eps") # This is Figure 3(b) par(mar=c(5,6,4,2)) rvec<-seq(from=0.05, to=1, by=0.05) c_star<-pmax(100*unlist(lapply(rvec, FUN=simplestati,b=0.25, v=.25/.67)),0) plot(c_star~rvec, type="l", main=c("Estimated c* (effect of drug on cancer probability)", "by level of loss aversion"), xlab="level of loss aversion (gamma)", ylab="", xlim=c(0, 1)) text(x=0.6, y=80, labels="act on null") text(x=0.23, y=40, labels="act on alternative") mtext(side=2, line=4, "minimum substantively meaningful relationship") mtext(side=2, line=3, "% change in cancer probability, drug - no drug (c*)") dev.copy2eps(file="drugrobust.eps") # This is Figure 3(a) ################################## # Replicate Figure 2 ################################## # draw a picture of a loss-averse utility function rm(list=ls()) setwd("C:/Dropbox/jesarey_documents/Substantive Robustness Project/Replication Materials") base<-function(a,x){(1/a)*(1-exp(-a*x));} x<-seq(from=-2, to=5, by=0.1) y<-base(0.3, seq(from=-2, to=5, by=0.1)) utility<-function(t,a,x){ k<-as.numeric((x-t)>0)*2-1 out<-a^(-k)*(x-t); return(out) } z<-utility(0,2.236068,x) par(mar=c(5,5,4,2)) plot(x,z,type="l", main="Loss Averse Utility Function", xlab=expression(beta), ylab=expression(paste("Utility from acting on alt. hypothesis (instead of null) given ", beta))) abline(a=0,b=1,lty=2) z<-utility(0,0.5,x) lines(x,z,lty=4, lwd=2) legend("bottomright", legend=c(expression(paste(gamma == 1, " ")),expression(gamma > 1), expression(gamma < 1)), lty=c(2,1,4), lwd=c(1,1,2)) dev.copy2eps(file="riskaversion.eps") # This is Figure 2