public

Visually-weighted regression plot, with Zelig

  • Download Gist
Visually_weighted_regression_Zelig.R
R
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50
# 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...

Please sign in to comment on this gist.

Something went wrong with that request. Please try again.