Skip to content

Embed URL

HTTPS clone URL

Subversion checkout URL

You can clone with
or
.
Download ZIP
Visually-weighted regression plot, with Zelig
# 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...
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Something went wrong with that request. Please try again.