#Overall setup rm(list = ls()) # delete all objects in memory library(foreign) # library for foreign (non-R) data files #library(arm) setwd("") # Define your working directory qra=read.dta("koenigmaeder_strategiccompliance_ajps.dta") attach(qra) #Model 1 #u11:znumber_of_provisions,zdelegationratio,zefficient_score #u13:const #u14:zmsdis_sum,zmss_sum,zcouconf_sum,typodir,zinterest,numissues #u23:0 #u24:const,zcomdis_sum,zcomss_sum,zcouconf_sum,typodir,zinterest,dummytime depvar=cbind(i,c,d) # dummies for dependent variable indvar=cbind( znumber_of_provisions,zdelegationratio,zefficient_score,zmsdis_sum,zmss_sum, zcouconf_sum,typodir,zinterest,numissues,zcomgain_sum, zcomss_sum,zcouconf_sum,typodir,zinterest,dummytime ) # log-likelihood procedure (NOT the negative of it!) llik=function(b,X,Y) { i=as.matrix(Y[,1]) # i, c, d are dep var dummies c=as.matrix(Y[,2]) d=as.matrix(Y[,3]) r=nrow(Y) o=rep(1,r) x11=as.matrix(cbind(X[,1],X[,2],X[,3])) x13=as.matrix(cbind(o)) x14=as.matrix(cbind(X[,4],X[,5],X[,6],X[,7],X[,8],X[,9])) x24=as.matrix(cbind(o,X[,10],X[,11],X[,12],X[,13],X[,14],X[,15])) b11=as.matrix(b[1:3]) b13=as.matrix(b[4]) b14=as.matrix(b[5:10]) b24=as.matrix(b[11:17]) u11=x11 %*% b11 u13=x13 %*% b13 u14=x14 %*% b14 u24=x24 %*% b24 p4=pnorm(u24/sqrt(2)) p3=1-p4 p2=pnorm((p3*u13+p4*u14-u11)/sqrt(p3^2+p4^2+1)) p1=1-p2 lnpi=log(p1) lnpc=log(p2*p3) lnpd=log(p2*p4) llik=i*lnpi+c*lnpc+d*lnpd sum(llik) } # MLE using optim() stval=runif(17)*0 strat.mle=optim(stval,llik,hessian=TRUE,method="BFGS",control=list(fnscale=-1,trace=1),X=indvar,Y=depvar) # Print results #u11:znumber_of_provisions,zdelegationratio,zefficient_score #u13:const #u14:zmsdis_sum,zmss_sum,zcouconf_sum,typodir,zinterest,numissues #u23:0 #u24:const,zcomdis_sum,zcomss_sum,zcouconf_sum,typodir,zinterest,dummytime Variable=c( "u11:Number of major provisions","u11:Delegation ratio","u11:Bureaucratic efficiency", "u13:Constant", "u14:Member state's disagreement","u14:Member State's saliency","u14:Diversity of member states' interests","u14:Type of directive (1 - CM & EP)","u14:Interest group pattern","u14:Number of issues", "u24:Constant","u24:Commission's agreement","u24:Commission's saliency","u24:Diversity of member states' interests","u24:Type of directive (1 - CM & EP)","u24:Interest group pattern","u24:Issue of timeliness") Parameter=as.matrix(strat.mle$par) StdError=as.matrix(sqrt(diag(solve(-1 * strat.mle$hessian)))) Tstat=as.matrix(Parameter/StdError) pvalue<-2*pt(abs(Tstat),length(i)-length(Parameter),lower.tail=FALSE) strat.mle.table=cbind(format(Variable,justify="left",width=20), prettyNum(Parameter,justify="right",digits=4), prettyNum(StdError,justify="right",digits=4), #prettyNum(Tstat,justify="right",digits=4), prettyNum(pvalue,justify="right",digits=4)) cat("\n") #Results=c("Variable","Parameter","StdError","t","pvalue"); Results=c("Variable","Parameter","StdError","pvalue"); print(format(rbind(Results,strat.mle.table),width=12,justify="right")) cat("\nLog-Likelihood = ",strat.mle$value,"\n") cat("N = ",length(i),"\n\n") # Predicted Probabilities n=299 const=rep(1,n) x11=cbind(znumber_of_provisions,zdelegationratio,zefficient_score) x13=const x14=cbind(zmsdis_sum,zmss_sum,zcouconf_sum,typodir,zinterest,numissues) x24=cbind(const,zcomgain_sum,zcomss_sum,zcouconf_sum,typodir,zinterest,dummytime) b11=matrix(strat.mle$par[1:3],3,1) b13=matrix(strat.mle$par[4],1,1) b14=matrix(strat.mle$par[5:10],6,1) b24=matrix(strat.mle$par[11:17],7,1) u11=x11 %*% b11 u13=x13 %*% b13 u14=x14 %*% b14 u24=x24 %*% b24 #Probabilities p4=pnorm(u24/sqrt(2)) p3=1-p4 p2=pnorm((p3*u13+p4*u14-u11)/sqrt(p3^2+p4^2+1)) p1=1-p2 pimpl=p1 pcap=p2*p3 pdis=p2*p4 #Create new variable depend pdepend <- rep(NA,299) pdepend[pcap>pimpl & pcap>pdis] <- 3 pdepend[pdis>pimpl & pdis>pcap] <- 4 pdepend[pimpl>pcap & pimpl>pdis] <- 1 #Comparison between observed and predicted crosstable <- table(pdepend, depend) crosstable #Corretly predicted n=299 cpredicted=(crosstable[1,1]+crosstable[2,2]+crosstable[3,3])/n print(cpredicted) #DIFFERENT OUTCOME PROBABILITIES n=299 const=rep(1,n) #Summary of independent variables summary(zmsdis_sum) summary(zmss_sum) summary(znumber_of_provisions) summary(zdelegationratio) summary(zefficient_score) summary(zcouconf_sum) summary(typodir) summary(zinterest) summary(numissues) summary(zcomgain_sum) summary(zcomss_sum) summary(dummytime) #Low values of indvar (1rd Qu.) low_zmsdis_sum=rep(-0.6616,n) low_zmss_sum=rep(-0.6626,n) low_znumber_of_provisions=rep(-0.6192,n) low_zdelegationratio=rep(-0.9686,n) low_zefficient_score=rep(-0.344,n) low_zcouconf_sum=rep(-0.7757,n) low_typodir=rep(1,n) low_zinterest=rep(-0.789,n) low_numissues=rep(1,n) low_zcomgain_sum=rep(-0.7757,n) low_zcomss_sum=rep(-0.7642,n) low_dummytime=rep(0,n) #Mean values of indvar mean_zmsdis_sum=rep(0,n) mean_zmss_sum=rep(0,n) mean_znumber_of_provisions=rep(0,n) mean_zdelegationratio=rep(0,n) mean_zefficient_score=rep(0,n) mean_zcouconf_sum=rep(0,n) mean_typodir=rep(1,n) mean_zinterest=rep(0,n) mean_numissues=rep(1.906,n) mean_zcomgain_sum=rep(0,n) mean_zcomss_sum=rep(0,n) mean_dummytime=rep(0,n) #Moderate values of indvar (3rd Qu.) mode_zmsdis_sum=rep(0.3854,n) mode_zmss_sum=rep(0.6348,n) mode_znumber_of_provisions=rep(0.1397,n) mode_zdelegationratio=rep(0.6642,n) mode_zefficient_score=rep(0.7768,n) mode_zcouconf_sum=rep(0.4183,n) mode_typodir=rep(1,n) mode_zinterest=rep(0.8686,n) mode_numissues=rep(2,n) mode_zcomgain_sum=rep(0.4586,n) mode_zcomss_sum=rep(0.7777,n) mode_dummytime=rep(0,n) #X-VARIABLE: DIVERSITY OF STATES' INTERESTS #zcouconf_sum n=299 const=rep(1,n) x=seq(-2,2.7,length=n) # Pr(impl/dis/cap|zcouconf_sum) holding other indvar at min/mean/max values #u11:znumber_of_provisions,zdelegationratio,zefficient_score #u13:const #u14:zmsdis_sum,zmss_sum,zcouconf_sum,typodir,zinterest,numissues #u23:0 #u24:const,zcomdis_sum,zcomss_sum,zcouconf_sum,typodir,zinterest,dummytime # Predicted Probabilities for LOW VALUES n=299 const=rep(1,n) x11=cbind(low_znumber_of_provisions,low_zdelegationratio,low_zefficient_score) x13=const x14=cbind(low_zmsdis_sum,low_zmss_sum,x,low_typodir,low_zinterest,low_numissues) x24=cbind(const,low_zcomgain_sum,low_zcomss_sum,x,low_typodir,low_zinterest,low_dummytime) b11=matrix(strat.mle$par[1:3],3,1) b13=matrix(strat.mle$par[4],1,1) b14=matrix(strat.mle$par[5:10],6,1) b24=matrix(strat.mle$par[11:17],7,1) u11=x11 %*% b11 u13=x13 %*% b13 u14=x14 %*% b14 u24=x24 %*% b24 #Probabilities p4=pnorm(u24/sqrt(2)) p3=1-p4 p2=pnorm((p3*u13+p4*u14-u11)/sqrt(p3^2+p4^2+1)) p1=1-p2 lowpimpl=p1 lowpcap=p2*p3 lowpdis=p2*p4 # Predicted Probabilities for MEAN VALUES n=299 const=rep(1,n) x11=cbind(mean_znumber_of_provisions,mean_zdelegationratio,mean_zefficient_score) x13=const x14=cbind(mean_zmsdis_sum,mean_zmss_sum,x,mean_typodir,mean_zinterest,mean_numissues) x24=cbind(const,mean_zcomgain_sum,mean_zcomss_sum,x,mean_typodir,mean_zinterest,mean_dummytime) b11=matrix(strat.mle$par[1:3],3,1) b13=matrix(strat.mle$par[4],1,1) b14=matrix(strat.mle$par[5:10],6,1) b24=matrix(strat.mle$par[11:17],7,1) u11=x11 %*% b11 u13=x13 %*% b13 u14=x14 %*% b14 u24=x24 %*% b24 #Probabilities p4=pnorm(u24/sqrt(2)) p3=1-p4 p2=pnorm((p3*u13+p4*u14-u11)/sqrt(p3^2+p4^2+1)) p1=1-p2 meanpimpl=p1 meanpcap=p2*p3 meanpdis=p2*p4 # Predicted Probabilities for MODERATE VALUES n=299 const=rep(1,n) x11=cbind(mode_znumber_of_provisions,mode_zdelegationratio,mode_zefficient_score) x13=const x14=cbind(mode_zmsdis_sum,mode_zmss_sum,x,mode_typodir,mode_zinterest,mode_numissues) x24=cbind(const,mode_zcomgain_sum,mode_zcomss_sum,x,mode_typodir,mode_zinterest,mode_dummytime) b11=matrix(strat.mle$par[1:3],3,1) b13=matrix(strat.mle$par[4],1,1) b14=matrix(strat.mle$par[5:10],6,1) b24=matrix(strat.mle$par[11:17],7,1) u11=x11 %*% b11 u13=x13 %*% b13 u14=x14 %*% b14 u24=x24 %*% b24 #Probabilities p4=pnorm(u24/sqrt(2)) p3=1-p4 p2=pnorm((p3*u13+p4*u14-u11)/sqrt(p3^2+p4^2+1)) p1=1-p2 modepimpl=p1 modepcap=p2*p3 modepdis=p2*p4 #Figure 4a #Plot Pr(enforced compliance|zcouconf_sum) holding other indvar at min/mean/max values plot(x,meanpdis,ylab="Pr(Enforced Compliance)",xlab="Diversity of states' interests",ylim=c(0,1.0), type="l", lty="solid",lwd="2.0",cex.lab=1.5) lines(x,lowpdis,lty="F8",lwd="2.0") lines(x,modepdis,lty="dotted",lwd="2.0") legend(-2.0, 1.0, c("Mean", "Low", "Moderate"), cex=1, lty=c("solid","F8","dotted"), lwd=1.9) #Figure 4b #Plot Pr(compliance deficit|zcouconf_sum) holding other indvar at min/mean/max values plot(x,meanpcap,ylab="Pr(Deficit)",xlab="Diversity of states' interests",ylim=c(0,1.0), type="l", lty="solid",lwd="2.0",cex.lab=1.5) lines(x,lowpcap,lty="F8",lwd="2.0") lines(x,modepcap,lty="dotted",lwd="2.0") legend(-2.0, 1.0, c("Mean", "Low", "Moderate"), cex=1.5, lty=c("solid","F8","dotted"), lwd=1.0) #X-VARIABLE: MEMBER STATES' DISAGREEMENT #zmsdis_sum n=299 const=rep(1,n) x=seq(-1.2,4.2,length=n) # Pr(impl/dis/cap|zcouconf_sum) holding other indvar at min/mean/max values #u11:znumber_of_provisions,zdelegationratio,zefficient_score #u13:const #u14:zmsdis_sum,zmss_sum,zcouconf_sum,typodir,zinterest,numissues #u23:0 #u24:const,zcomdis_sum,zcomss_sum,zcouconf_sum,typodir,zinterest,dummytime # Predicted Probabilities for LOW VALUES n=299 const=rep(1,n) x11=cbind(low_znumber_of_provisions,low_zdelegationratio,low_zefficient_score) x13=const x14=cbind(x,low_zmss_sum,low_zcouconf_sum,low_typodir,low_zinterest,low_numissues) x24=cbind(const,low_zcomgain_sum,low_zcomss_sum,low_zcouconf_sum,low_typodir,low_zinterest,low_dummytime) b11=matrix(strat.mle$par[1:3],3,1) b13=matrix(strat.mle$par[4],1,1) b14=matrix(strat.mle$par[5:10],6,1) b24=matrix(strat.mle$par[11:17],7,1) u11=x11 %*% b11 u13=x13 %*% b13 u14=x14 %*% b14 u24=x24 %*% b24 #Probabilities p4=pnorm(u24/sqrt(2)) p3=1-p4 p2=pnorm((p3*u13+p4*u14-u11)/sqrt(p3^2+p4^2+1)) p1=1-p2 lowpimpl=p1 lowpcap=p2*p3 lowpdis=p2*p4 # Predicted Probabilities for MEAN VALUES n=299 const=rep(1,n) x11=cbind(mean_znumber_of_provisions,mean_zdelegationratio,mean_zefficient_score) x13=const x14=cbind(x,mean_zmss_sum,mean_zcouconf_sum,mean_typodir,mean_zinterest,mean_numissues) x24=cbind(const,mean_zcomgain_sum,mean_zcomss_sum,mean_zcouconf_sum,mean_typodir,mean_zinterest,mean_dummytime) b11=matrix(strat.mle$par[1:3],3,1) b13=matrix(strat.mle$par[4],1,1) b14=matrix(strat.mle$par[5:10],6,1) b24=matrix(strat.mle$par[11:17],7,1) u11=x11 %*% b11 u13=x13 %*% b13 u14=x14 %*% b14 u24=x24 %*% b24 #Probabilities p4=pnorm(u24/sqrt(2)) p3=1-p4 p2=pnorm((p3*u13+p4*u14-u11)/sqrt(p3^2+p4^2+1)) p1=1-p2 meanpimpl=p1 meanpcap=p2*p3 meanpdis=p2*p4 # Predicted Probabilities for MODERATE VALUES n=299 const=rep(1,n) x11=cbind(mode_znumber_of_provisions,mode_zdelegationratio,mode_zefficient_score) x13=const x14=cbind(x,mode_zmss_sum,mode_zcouconf_sum,mode_typodir,mode_zinterest,mode_numissues) x24=cbind(const,mode_zcomgain_sum,mode_zcomss_sum,mode_zcouconf_sum,mode_typodir,mode_zinterest,mode_dummytime) b11=matrix(strat.mle$par[1:3],3,1) b13=matrix(strat.mle$par[4],1,1) b14=matrix(strat.mle$par[5:10],6,1) b24=matrix(strat.mle$par[11:17],7,1) u11=x11 %*% b11 u13=x13 %*% b13 u14=x14 %*% b14 u24=x24 %*% b24 #Probabilities p4=pnorm(u24/sqrt(2)) p3=1-p4 p2=pnorm((p3*u13+p4*u14-u11)/sqrt(p3^2+p4^2+1)) p1=1-p2 modepimpl=p1 modepcap=p2*p3 modepdis=p2*p4 #Figure 4c #Plot Pr(dis|zmsdis_sum) holding other indvar at min/mean/max values plot(x,meanpdis,ylab="Pr(Enforced Compliance)",xlab="Member state's disagreement",ylim=c(0,1.0), type="l", lty="solid",lwd="2.0",cex.lab=1.5) lines(x,lowpdis,lty="F8",lwd="2.0") lines(x,modepdis,lty="dotted",lwd="2.0") legend(-1.2, 1.0, c("Mean", "Low", "Moderate"), cex=1.2, lty=c("solid","F8","dotted"), lwd=1.0) #Figure 4d #Plot Pr(cap|zmsdis_sum) holding other indvar at min/mean/max values plot(x,meanpcap,ylab="Pr(Deficit)",xlab="Member state's disagreement",ylim=c(0,1.0), type="l", lty="solid",lwd="2.0",cex.lab=1.5) lines(x,lowpcap,lty="F8",lwd="2.0") lines(x,modepcap,lty="dotted",lwd="2.0") legend(-1.2, 1.0, c("Mean", "Low", "Moderate"), cex=1.2, lty=c("solid","F8","dotted"), lwd=1.0) #Figure 5 #3D Plots outcome probabilies for a change in two explanatory variables #Outcome function #x: Diversity of states range(zcouconf_sum) x <- seq(-1.347835,2.570057,length=299) #y: MS disagreement range(zmsdis_sum) y <- seq(-0.9607449,4.0349951,length=299) outcomedis <- function(x,y){ n=299 const=rep(1,n) x11=cbind(mean_znumber_of_provisions,mean_zdelegationratio,mean_zefficient_score) x13=const x14=cbind(y,mean_zmss_sum,x,mean_typodir,mean_zinterest,mean_numissues) x24=cbind(const,mean_zcomgain_sum,mean_zcomss_sum,x,mean_typodir,mean_zinterest,mean_dummytime) b11=matrix(strat.mle$par[1:3],3,1) b13=matrix(strat.mle$par[4],1,1) b14=matrix(strat.mle$par[5:10],6,1) b24=matrix(strat.mle$par[11:17],7,1) u11=x11 %*% b11 u11 <- as.numeric(u11) u13=x13 %*% b13 u13 <- as.numeric(u13) u14=x14 %*% b14 u14 <- as.numeric(u14) u24=x24 %*% b24 u24 <- as.numeric(u24) #Probabilities p4=pnorm(u24/sqrt(2)) p3=1-p4 p2=pnorm((p3*u13+p4*u14-u11)/sqrt(p3^2+p4^2+1)) p1=1-p2 pimpl=p1 pcap=p2*p3 pdis=p2*p4 return(pdis) } zdis <- outer(x,y,outcomedis) persp(x, y, zdis, theta = 30, phi = 30, expand = 0.9, col = "gray80", border= NA, shade = 0.8, xlab="Diversity of states' interests",ylab="MS disagreement",zlab="Pr(Enforced Compliance)", ticktype="detailed",cex.lab=1,zlim=c(0,1))