#################### # Create relogit predicted probabilities using Zelig and ggplot2 # Two Sword Lengths: Losers' Consent and Violence in National Legislatures (Working Paper 2012) # Christopher Gandrud # Updated 26 April 2012 ################### ## Load required packages library(RCurl) library(Zelig) library(ggplot2) library(reshape2) library(gridExtra) ## Load data leg <- read.csv("http://dl.dropbox.com/u/12581470/code/Replicability_code/Leg_Violence/Gandrud_leg_violence_ggplot_example.csv") ## List of variables included in the rare event logistic regression vars <- c("country", "year", "violence", "DemAge", "maj", "HighProp", "tenshort") ## Subset data to only include complete cases leg.comp <- leg[complete.cases(leg[vars]),] ## Run rare event logistic regression. ## Note the regression uses prior correction (tau) with the full sample and 72 incidences of observed violence (see King and Zeng 2002). Also I used WEAVE robust standard errors (see Lumley and Heagrty 1999). Model <- zelig(violence ~ DemAge + maj + HighProp + tenshort, model = "relogit", data = leg.comp, tau = 72/4224, robust = list(method = "weave")) ## Ranges of fitted values HighProp.r <- c(0, 1) dem.r <- 0:85 maj.r <- 20:100 ############# Create Individual Variable Plots ################# ### Age of Democracy ### # Set fitted values Model.DemAge <-setx(Model, DemAge = dem.r) # Simulate quantities of interest Model.DemSim <- sim(Model, x = Model.DemAge) # Extract expected values from simulations Model.demAge.e <- Model.DemSim$qi Model.demAge.e <- data.frame(Model.demAge.e$ev) Model.demAge.e <- melt(Model.demAge.e, measure = 1:86) # Remove "X" from variable Model.demAge.e$variable <- as.numeric(gsub("X", "", Model.demAge.e$variable)) # Plot Model.demAge.p <- ggplot(Model.demAge.e, aes(variable, value)) + geom_point(shape = 21, color = "gray30", alpha = I(0.05)) + stat_smooth() + scale_x_continuous(breaks = c(0, 21, 51, 85), labels = c("0", "20", "50", "85")) + scale_y_continuous(breaks = c(0, 0.02, 0.05, 0.1, 0.15), limits = c(0, 0.12)) + xlab("\nAge of Democracy") + ylab("") + opts(axis.title.x = theme_text(size = 12, vjust = 0)) + theme_bw() ### High Proportionality ### ## Set fitted values Model.HighProp <- setx(Model, HighProp = HighProp.r) ## Simulate quantities of interest Model.HighPropSim <- sim(Model, x = Model.HighProp) ## Extract expected values from simulations Model.HighProp.e <- (Model.HighPropSim$qi) Model.HighProp.e <- data.frame(Model.HighProp.e$ev) Model.HighProp.e <- melt(Model.HighProp.e, measure = 1:2) # Remove "X" from variable Model.HighProp.e$variable <- as.numeric(gsub("X", "", Model.HighProp.e$variable)) # Plot Model.HighProp.p <- ggplot(Model.HighProp.e, aes(variable, value)) + geom_point(shape = 21, color = "gray30", alpha = I(0.05)) + stat_smooth(method = "lm", se = FALSE) + scale_x_reverse(breaks = c(1, 2), labels = c("Higher", "Very Low")) + scale_y_continuous(breaks = c(0, 0.02, 0.05, 0.1, 0.15), labels = c("", "", "", "", ""), limits = c(0, 0.12)) + xlab("\nDisproportionality") + ylab("") + opts(axis.title.x = theme_text(size = 12, vjust = 0)) + theme_bw() ### Majority ### ## Set fitted values Model.maj1 <-setx(Model, maj = maj.r) ## Simulate quantities of interest Model.majSim <- sim(Model, x = Model.maj1) ## Extract expected values from simulations Model.maj.e <- (Model.majSim$qi) Model.maj.e <- data.frame(Model.maj.e$ev) Model.maj.e <- melt(Model.maj.e, measure = 1:81) ## Remove "X" from variable Model.maj.e$variable <- as.numeric(gsub("X", "", Model.maj.e$variable)) ## Put in terms of the original variable percentage Model.maj.e$variable = Model.maj.e$variable + 19 ## Plot Model.maj.p <- ggplot(Model.maj.e, aes(variable, value)) + geom_point(shape = 21, color = "gray30", alpha = I(0.05)) + stat_smooth(method = "lm", se = FALSE) + scale_y_continuous(breaks = c(0, 0.02, 0.05, 0.1, 0.15), labels = c("", "", "", "", ""), limits = c(0, 0.12)) + xlab("\nGovernment Majority") + ylab("") + theme_bw() ############# Combine Individual Variable Plots ################# predicted.combine <- grid.arrange(Model.demAge.p, Model.HighProp.p, Model.maj.p, ncol = 3, left = "Predicted Probability of Violence") print(predicted.combine)