Visually-weighted regression plot, with Zelig

Visually_weighted_regression_Zelig.R

# A simple approach to visually-weighted regression plots, with Zelig

doInstall <- TRUE  # Change to FALSE if you don't want packages installed.
toInstall <- c("ggplot2", "reshape2", "Zelig")
if(doInstall){install.packages(toInstall, repos = "http://cran.us.r-project.org")}
lapply(toInstall, library, character.only = TRUE)

# Generate some data:
nn <- 1000
myData <- data.frame(X = rnorm(nn),
                     Binary = sample(c(0, 1), nn, replace = T))
myData$Y <- with(myData, X * 2 + Binary * 10 + rnorm(nn, sd = 10))
myData$Y <- (myData$Y > 1) * 1
head(myData)

# Model, using Zelig
myModel <- zelig(Y ~ X + Binary, model="logit", data=myData)

# Simulate from model, using Zelig
rangeX <- with(myData, seq(min(X, na.rm = T), max(X, na.rm = T), length.out = 100))
lowScenario <- setx(myModel, X = rangeX, Binary = 0)
highScenario <- setx(myModel, X = rangeX, Binary = 1)
simulatedScenarios <- sim(myModel, x = lowScenario, x1 = highScenario)

# Make a tall data.frame from the Zelig simulation object
toReshape <- data.frame(x = simulatedScenarios[[2]][, "X"],
                        t(simulatedScenarios[[6]]$ev))
longEV <- melt(toReshape, id.vars = "x")

toReshape2 <- data.frame(x = simulatedScenarios[[2]][, "X"],
                         t(simulatedScenarios[[6]]$fd) + t(simulatedScenarios[[6]]$ev))
longFD <- melt(toReshape2, id.vars = "x")

longEV$Setting <- "Minimum"
longFD$Setting <- "Maximum"

longEV <- data.frame(rbind(longEV, longFD))

# Plot
zp1 <- ggplot(data = longEV,
              aes(x = x, y = value, group = paste(variable, Setting), colour = factor(Setting)))
zp1 <- zp1 + geom_line(alpha = I(1/sqrt(nrow(simulatedScenarios[[6]]$ev))))
zp1 <- zp1 + scale_x_continuous("X-axis label", expand = c(0, 0))
zp1 <- zp1 + scale_y_continuous("Y-axis label", limits = c(0, 1), expand = c(0, 0))
zp1 <- zp1 + scale_colour_brewer(palette="Set1", labels = c("High", "Low"))  # Change colour palette
zp1 <- zp1 + guides(colour = guide_legend("First-difference\nvariable",
                                          override.aes = list(alpha = 1)))  # Avoid an alpha-related legend problem
zp1 <- zp1 + ggtitle("Plot title")
zp1 <- zp1 + theme_bw()
print(zp1)  # This might take a few seconds...