Code for post on physical integrity rights and democracy

Malevolent democracies and benevolent autocracies

library(ggplot2)
library(plyr)
library(maps)

CIRI <- read.csv("../Data/CIRI Human Rights Data.csv")
codes <- read.csv("../Data/extended codes.csv")
polity <- read.csv("../Data/p4v2011.csv")
#uds <- read.csv("../Data/uds_summary.csv")
#utip <- read.csv("../Data/EHII2008.csv")
pwt <- read.csv("../Data/pwt70_w_country_names.csv")
pts <- read.csv("../Data/PTS_2011.csv")

extraYears <- rbind(codes[ codes$year ==2008, ],codes[ codes$year ==2008, ], codes[ codes$year ==2008, ])
extraYears$year <- rep(c(2009,2010,2011),each=190)
codes <- rbind(codes,extraYears) 
rm(extraYears)

names(pts)[3] <- "ccode"
names(pts)[5] <- "year"

pts[ pts$Country == "USSR", ]$ccode <- 364
ptspluspolity <- merge(pts,polity,by=c("ccode","year"))
ptspluspolity$predicted <- -0.2*ptspluspolity$polity2 + 3
ptspluspolity$difference <- abs(ptspluspolity$avg - ptspluspolity$predicted)
ptspluspolity$ccode <- gsub("364","365",ptspluspolity$ccode,fixed=TRUE)

ptspluspolity <- merge(ptspluspolity,codes, by=c("ccode","year"))
ptspluspolity$polityCut <- polityCut <- cut(ptspluspolity$polity,breaks=c(-100,-11,-6,6,10),labels=c("Interruption","Autocracy","Anocracy","Democracy"))
pts[ pts$Country == "USSR", ]$ccode <- 365

names(CIRI)[2] <- "year"
CIRI$ccode <- CIRI$COW


CIRI.2 <- merge(CIRI,codes,by=c("year","ccode"))


names(polity)[2] <- "politycode"
names(CIRI.2)[6] <- "politycode"

CIRIplusPolity <- merge(CIRI.2,polity,by=c("politycode","year"))


CIRIplusPolity$polityCut <- cut(CIRIplusPolity$polity,breaks=c(-100,-11,-6,6,10),labels=c("Interruption","Autocracy","Anocracy","Democracy"))

CIRIplusPolity$predicted <- CIRIplusPolity$polity2 * 0.4 + 4
CIRIplusPolity$difference <- abs(CIRIplusPolity$PHYSINT - CIRIplusPolity$predicted)
CIRIplusPolity$polity2.positive <- CIRIplusPolity$polity2 + 10
CIRIplusPolity$difference2 <- CIRIplusPolity$PHYSINT - CIRIplusPolity$polity2.positive * (8/20)
CIRIplusPolity$benevolenceCat <- cut(CIRIplusPolity$difference2,breaks=c(-9,-6,-3,3,6,8),labels=c("Highly malevolent (< -6)","Mildly malevolent (< -3)","Neutral (< 3)","Mildly benevolent (< 6)","Highly benevolent (<= 8)"))


names(CIRI)[5] <- "politycode"
dd <- read.csv("../Data/ddrevisited-data-v1.csv")
CIRIplusDD <- merge(CIRI,dd, by=c("politycode","year"))
CIRIPolityDD <- merge(CIRIplusPolity,CIRIplusDD, by=c("politycode","year"))


#Figure 1: Global Mean of the CIRI Index of Physical Integrity Rights
fig.1 <- qplot(data=CIRIplusPolity, x= year, y= PHYSINT,stat="summary",fun.y=mean)+ylim(0,8)+geom_smooth()+labs(y="Global Mean of CIRI Physical Integrity Rights Index")
fig.1
ggsave(plot=fig.1,"Figure 1 Global Mean of Physical Integrity Rights.png")

#Figure 2: Global Mean of Political Terror Scale
fig.2 <- qplot(data=ptspluspolity, x= year, y= avg,stat="summary",fun.y=mean)+ylim(1,5)+geom_smooth()+labs(y="Global Mean of Political Terror Scale")
fig.2
ggsave(plot=fig.2,"Figure 2 Global Mean of Political Terror Scale.png")

#Figure 3: Global Mean of Polity2 Score
fig.3 <- qplot(data=ptspluspolity, x= year, y= polity2,stat="summary",fun.y=mean)+ylim(-10,10)+geom_smooth()+labs(y="Global Mean of Polity 2 Score")
fig.3
ggsave(plot=fig.3,"Figure 3 Global Mean of Polity 2 Score.png")

#Figure 4: Global Trends by regime
fig.4 <- qplot(data=CIRIplusPolity,x=year,y=PHYSINT,geom="point",position="jitter",colour=polityCut)+geom_smooth()+labs(y="Physical integrity rights index",colour="Polity category")
fig.4
ggsave(plot=fig.4,"Figure 4 Global trends in PHYSINT by regime.png")

#Figure 5: Physical Integrity Index, by regime category
fig.5 <- qplot(data=CIRIplusPolity[ !is.na(CIRIplusPolity$PHYSINT) & !is.na(CIRIplusPolity$polityCut) , ],x=PHYSINT,geom="histogram",fill=polityCut,position="fill")+labs(x="CIRI Index of Physical Integrity Rights",fill="Polity category", y="Proportion of states")
fig.5
ggsave(plot=fig.5,"Figure 5 Physical integrity by regime category.png")

#Figure 6: Physical Integrity Index, by regime category (DD data)
fig.6 <- qplot(data=CIRIplusDD,x=PHYSINT,geom="histogram",fill=factor(democracy, labels=c("Dictatorship","Democracy")),position="fill")+labs(x="CIRI Index of Physical Integrity Rights",fill="DD regime", y="Proportion of states")
fig.6
ggsave(plot=fig.6,"Figure 6 Physical integrity by regime category.png")

#Figure 7: Physical Integrity Index, by regime category (DD data)
fig.7 <- qplot(data=CIRIplusDD,x=factor(regime, labels=c("Parliamentary democracy","Mixed Pres/Parl democracy","Presidential democracy","Civilian Dictatorship","Military Dictatorship","Monarchy")),geom="histogram",fill=as.factor(PHYSINT))+labs(x="DD regime category",fill="Phys. Int. Index", y="Number of country-years in dataset")+scale_fill_brewer(type="div",palette="RdBu",na.value="grey")+coord_flip()
fig.7
ggsave(plot=fig.7,"Figure 7 Physical integrity by regime category.png")

#Figure 8: Geographical distribution of regimes
worldMapData <- map("world", plot=FALSE, fill=TRUE)
worldMapData <- fortify(worldMapData)
worldMap <- qplot(data=worldMapData,x=long,y=lat,geom="path",group=group)

CIRIplusPolity.2 <- subset(CIRIplusPolity, !is.na(PHYSINT))
CIRIplusPolity.2 <- ddply(CIRIplusPolity.2, .(politycode), transform, yearsMeasured = length(PHYSINT))
CIRIplusPolity.2 <- ddply(CIRIplusPolity.2, .(politycode), transform, meanPHYSINT = mean(PHYSINT))
CIRIplusPolity.2 <- ddply(CIRIplusPolity.2, .(politycode), transform, meanPolity2 = mean(polity2,na.rm=TRUE))

CIRIplusPolity.2 <- unique(CIRIplusPolity.2[ , c("politycode","CTRY","CAPLONG","CAPLAT","dpicode","yearsMeasured","meanPHYSINT","meanPolity2")])
CIRIplusPolity.2[ CIRIplusPolity.2$yearsMeasured > 30, ]$yearsMeasured <- 30
fig.8 <- worldMap + geom_point(data=CIRIplusPolity.2,
                           aes(x=CAPLONG, 
                               y= CAPLAT,
                               colour = meanPHYSINT,
                               size=yearsMeasured,
                               group=NA,
                               shape=cut(0-meanPolity2,breaks=c(-11,-6,5,10),labels=c("(6,10] - Democracy", "(-6,6] - Anocracy","[-10,-6] - Autocracy"))))

fig.8 <- fig.8 + scale_colour_gradient2(low="red",high="blue",midpoint=4)

fig.8 <- fig.8 + labs(colour = "Physical Integrity Index",size = "Years with data", shape="Avg. polity2 category",title="Benevolent Autocracies and Malevolent Democracies, 1981-2010")

#Label outliers
benevolentAutocracies <- CIRIplusPolity.2[ CIRIplusPolity.2$meanPolity2 <= 6 & CIRIplusPolity.2$meanPHYSINT >= 6, ]
malevolentDemocracies <- CIRIplusPolity.2[ CIRIplusPolity.2$meanPolity2 > 6 & CIRIplusPolity.2$meanPHYSINT < 4, ]

fig.8 <- fig.8 + geom_text(data=benevolentAutocracies,aes(x=CAPLONG, y= CAPLAT,group=NA, fontface=3,label=benevolentAutocracies$dpicode))

fig.8 <- fig.8 + geom_text(data=malevolentDemocracies,aes(x=CAPLONG, y= CAPLAT,group=NA, fontface=1,label=malevolentDemocracies$dpicode))
fig.8
ggsave(plot=fig.8,file="Figure 8 Benevolent Autocracies and Malevolent Democracies 1981-2010.png",width=17, height=10)

# Fig. 9: Graph of Ciri for USA
fig.9 <- qplot(data=CIRIplusPolity[ CIRIplusPolity$ccode == 2, ], x = year, y = PHYSINT, size=polity2)+ylim(0,8)+ geom_smooth() + labs(y="Physical Integrity Index",size="Polity2 score",title="Physical Integrity Index in the USA, 1981-2010")+geom_vline(xintercept=2001,colour="red")
ggsave(plot=fig.9,file="Figure 9 Physical Integrity Index in the USA.png")

#Figure 10: distribution of countries by benevolence

cpolity <- subset(CIRIplusPolity, select=c("CIRI","year","polity2","polity","exconst","CAPLONG","CAPLAT"))
cddpolity <- merge(CIRIplusDD,cpolity)
rm(cpolity)

cddpolity <- subset(cddpolity, !is.na(PHYSINT) & exconst > 0)
cddpolity$benevolence <- sqrt((8 - cddpolity$exconst)*(cddpolity$PHYSINT + 1))
cddpolity$democracy2 <- factor(cddpolity$democracy,labels=c("Non-democracy","Democracy"))
cddpolity$polityCut <- cut(cddpolity$polity2, breaks=c(-11,-6,6,10),labels=c("Autocracy","Anocracy","Democracy"))
cddpolity$regime2 <- factor(cddpolity$regime,labels=c("Parl. Dem.","Mixed Dem.","Pres. Dem.","Civilian Dict.","Military Dict.","Monarchy"))

cddpolity <- ddply(cddpolity, .(CIRI), transform, yearsMeasured = length(PHYSINT))

cddpolity <- ddply(cddpolity, .(CIRI, regime2), transform, yearsMeasuredReg = length(PHYSINT))

cddpolity <- ddply(cddpolity, .(CIRI, democracy2), transform, yearsMeasuredDem = length(PHYSINT))

cddpolity <- ddply(cddpolity, .(CIRI, regime2), transform, meanBenevolenceReg = mean(benevolence,na.rm=TRUE))

cddpolity <- ddply(cddpolity, .(CIRI, democracy2), transform, meanBenevolenceDem = mean(benevolence,na.rm=TRUE))

cddpolity <- unique(cddpolity[ , c("CIRI","CTRY","CAPLONG","CAPLAT","dpicode","yearsMeasuredReg","yearsMeasured","meanBenevolenceReg","meanBenevolenceDem","regime2","democracy2")])
cddpolity[ cddpolity$yearsMeasured > 30, ]$yearsMeasured <- 28
cddpolity[ cddpolity$yearsMeasuredReg > 30, ]$yearsMeasuredReg <- 28


fig.10 <- qplot(data=cddpolity,x=reorder(CTRY,meanBenevolenceReg),y=meanBenevolenceDem,shape=regime2, size=yearsMeasuredReg,colour=regime2)+coord_flip()+facet_wrap(~democracy2)+geom_hline(yintercept=4,color="red")
fig.10 <- fig.10+ labs(y="Avg. Degree of Benevolence, 1981-2008", x="Country", shape="Regime type",colour="Regime type",size="Years \nwith \ndata")
fig.10
ggsave(plot=fig.10,file="Figure 10 Regimes by level of benevolence.png",height=25)
#Slideshow 1

library(cshapes)

rights.deviation.map <- function(year, df=CIRIplusPolity) {
  #   If using cshapes uncomment below
  
    if(year > 2008) {
      date <- as.Date("2008-1-1")
    }
    else {
      date <- as.Date(paste(year,"-1-1",sep=""))
    }
    
    worldMapData <- cshp(date=date)
    worldMapData <- fortify(worldMapData)
    worldMap <- qplot(data=worldMapData,x=long,y=lat,geom="path",group=group)
  
  data.year <- df[df$year == year, ]
  
  m <- worldMap + geom_point(data=data.year,
                             aes(x=CAPLONG, 
                                 y= CAPLAT,
                                 colour = difference2,
                                 group=NA,
                                 size=PHYSINT,
                                 shape=cut(0-polity,breaks=c(-11,-6,5,10,67,78,89),labels=c("Democracy","Anocracy","Autocracy","Interruption","Interregnum","Transition"))))
  
  m <- m + scale_colour_gradient2(low="red",high="blue",midpoint=0, labels=c("-8 - highly malevolent","-4","0","4","8 - highly benevolent"),limits=c(-8.01,8))
  
  m <- m + labs(colour = "Deviation from expectation", shape="Regime type",title=paste("State malevolence and state benevolence in the world, ",year), size="Phys. Integ. Index")
  m <- m + scale_size(range=c(3,8))
  
  #Label outliers
  
  m <- m + geom_text(data=data.year[ abs(data.year$difference2) > 5.5, ],aes(x=CAPLONG, y= CAPLAT,group=NA),label=data.year$dpicode[ abs(data.year$difference2) > 5.5 ],size=4)
  
  m
  
}

for(i in 1981:2010) {
  m <- rights.deviation.map(i)
  ggsave(m,file=paste("State malevolence and state benevolence in the World",i,".png"),width=19,height=10)
}


# Figure 11: the relationship between polity2 and PHYSINT over time
cp <- subset(CIRIplusPolity, !is.na(polity2) & !is.na(PHYSINT))
mods.cp <- function(df) {
  m <- lm(data=df, PHYSINT ~ polity2)
  m
} 
models.year <- dlply(cp, .(year), .fun=mods.cp)
coefs.year <- lapply(models.year, coef)
confints.year <- lapply(models.year, confint)
coefs.year <- ldply(coefs.year)
confints.year.polity2 <- ldply(confints.year)[ seq.int(from=2, by= 2, to=60), ]
names(coefs.year)[1] <- "year"
names(confints.year.polity2) <- c("year","low2.5","high97.5")
data <- merge(coefs.year, confints.year.polity2,by="year")
data$year <- as.numeric(data$year)
data$coldwar <- ifelse(data$year < 1990,TRUE,FALSE)
fig.11 <- qplot(data=data,x=year,y=polity2,group=coldwar)+geom_smooth(aes(ymin=low2.5,ymax=high97.5),data=data,stat="identity",colour=NA)+geom_vline(xintercept=1990,colour="red")+geom_smooth(method="lm",se=FALSE)+labs(y="Estimated strength of relationship between polity2")
fig.11
ggsave(fig.11,file="Figure 11 changes in the relationship.png")

#Figure 12: relationship between polity2 and PHYSINT, per country

models.ctry <- dlply(cp, .(CTRY), .fun=mods.cp)
coefs.ctry <- lapply(models.ctry, coef)
confints.ctry <- lapply(models.ctry, confint)
coefs.ctry <- ldply(coefs.ctry)
confints.ctry.polity2 <- ldply(confints.ctry)[ seq.int(from=2, by= 2, to=336), ]
names(coefs.ctry)[1] <- "CTRY"
names(confints.ctry.polity2) <- c("CTRY","low2.5","high97.5")
data <- merge(coefs.ctry, confints.ctry.polity2,by="CTRY")
cp <- ddply(cp, .(CTRY), transform, meanPHYSINT= mean(PHYSINT))
cp <- ddply(cp, .(CTRY), transform, meanPolity2= mean(polity2))
cp <- ddply(cp, .(CTRY), transform, meanPolity2Cut= cut(0-meanPolity2,breaks=c(-11,-5,5,10),labels=c("Democracy","Anocracy","Autocracy")))
cp <- ddply(cp, .(CTRY), transform, rangePolityChange= max(polity2)-min(polity2))
cp_extras <- unique(subset(cp, select=c("CTRY","un_continent_name","un_region_name","meanPHYSINT","meanPolity2","meanPolity2Cut","rangePolityChange")))
data <- merge(data,cp_extras,by="CTRY")
fig.12 <- qplot(data=data[complete.cases(data),],x=reorder(CTRY,polity2),y=polity2,ymax=high97.5,ymin=low2.5,colour=meanPHYSINT,shape=meanPolity2Cut,geom="pointrange")+coord_flip()+scale_colour_gradient2(high="blue",low="red",midpoint=4,limits=c(0,8))

fig.12 <- fig.12 + labs(x="Country",y="Strength of relationship between changes in polity2 score and changes in CIRI Index",colour="Mean \nPhysical \nIntegrity \nIndex",shape="Mean \nPolity2 \nCategory")+geom_hline(yintercept=0,colour="red")
fig.12
ggsave(plot=fig.12,"Figure 12 magnitude of relation between changes in polity2 and changes in CIRI per country.png",height=19)

# Figure 13
CIRIplusDD$regime2 <- factor(CIRIplusDD$regime,labels=c("Parl. Dem.","Mixed Dem.","Pres. Dem.","Civilian Dict.","Military Dict.","Monarchy"))

pwt <- read.csv("../Data/pwt70_w_country_names.csv")
library(countrycode)
pwt$isocode <-gsub("ROM","ROU",pwt$isocode)
pwt$isocode <-gsub("ZAR","COD",pwt$isocode)
pwt$isocode <-gsub("GER","DEU",pwt$isocode)

pwt$cowcode <- countrycode(pwt$isocode,"iso3c","cown")
pwt <- ddply(pwt,.(isocode), transform, cgdpgrowthch = 100*diff(rgdpch)/rgdpch)
pwt <- ddply(pwt,.(isocode), transform, cgdpgrowthl = 100*diff(rgdpl)/rgdpl)
cddpwt <- merge(CIRIplusDD,pwt,by=c("cowcode","year"))

fig.13 <- qplot(data=cddpwt,x=rgdpl,y=PHYSINT,log="x",colour=regime2,position="jitter")+geom_smooth(method="lm")
fig.13 <- fig.13 + labs(x="GDP per capita (log scale)", y="Physical Integrity Index",colour="Regime type")
fig.13
ggsave(plot=fig.13,"Figure 13 relationship between PHYSINT and gdp per capita by regime type.png", width=10)

# Fig. 14
library(reshape2)
utip <- read.csv("../Data/EHII2008.csv")
utip.melted <- melt(utip,id.vars=1:2)
utip.melted$year <- as.numeric(gsub("X","",utip.melted$variable))
utip.melted <- utip.melted[ !is.na(utip.melted$year), ]
names(utip.melted)[4] <- "gini"
names(utip.melted)[1] <- "country.utip"
names(utip.melted)[2] <- "dpicode"
utip.melted <- utip.melted[ , c(1:2,4:5)]
utip <- utip.melted
rm(utip.melted)

CIRIplusDD$democracy2 <- factor(CIRIplusDD$democracy,labels=c("Non-democracy","Democracy"))
cddutip <- merge(CIRIplusDD,utip,by=c("dpicode","year"))

fig.14 <- qplot(data=cddutip,x=gini,y=PHYSINT,colour=democracy2,position="jitter")+geom_smooth(method="lm")
fig.14 <- fig.14 + labs(x="GINI index of inequality", y="Physical Integrity Index",colour="Regime type")
fig.14
ggsave(plot=fig.14,"Figure 14 relationship between PHYSINT and inequality by regime type.png", width=10)

#Figure 15 GDP groth and repressiveness
fig.15 <- qplot(data=cddpwt,x=cgdpgrowthl,y=PHYSINT,colour=democracy2,position="jitter")+geom_smooth(method="lm")+xlim(-10,10)
fig.15 <- fig.15 + labs(x="Growth in GDP per capita (%)", y="Physical Integrity Index",colour="Regime type")
fig.15
ggsave(plot=fig.15,"Figure 15 relationship between PHYSINT and GDP growth by regime type.png", width=10)