Predict probability graphs with zelig and ggplot2

leg_violence_predict.R

####################
# 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)