Skip to content
Create a gist now

Instantly share code, notes, and snippets.

A coefficient plot for multiple models
doInstall <- TRUE
toInstall <- c("ggplot2")
if(doInstall){install.packages(toInstall, repos = "http://cran.us.r-project.org")}
lapply(toInstall, library, character.only = TRUE)
ANES <- read.csv("http://www.oberlin.edu/faculty/cdesante/assets/downloads/ANES.csv")
ANES <- ANES[ANES$year == 2008, -c(1, 11, 17)] # Limit to just 2008 respondents,
head(ANES) # remove some non-helpful variables
# Fit several models with the same DV:
model1 <- lm(pid7 ~ ideo7 + female + age + south, data = ANES)
model2 <- lm(pid7 ~ ideo7 + female + age + female:age, data = ANES)
model3 <- lm(pid7 ~ ideo7, data = ANES) # These are just arbitrary examples
# Put model estimates into temporary data.frames:
model1Frame <- data.frame(Variable = rownames(summary(model1)$coef),
Coefficient = summary(model1)$coef[, 1],
SE = summary(model1)$coef[, 2],
modelName = "South Indicator")
model2Frame <- data.frame(Variable = rownames(summary(model2)$coef),
Coefficient = summary(model2)$coef[, 1],
SE = summary(model2)$coef[, 2],
modelName = "Age Interaction")
model3Frame <- data.frame(Variable = rownames(summary(model3)$coef),
Coefficient = summary(model3)$coef[, 1],
SE = summary(model3)$coef[, 2],
modelName = "Univariate")
# Combine these data.frames
allModelFrame <- data.frame(rbind(model1Frame, model2Frame, model3Frame)) # etc.
# Specify the width of your confidence intervals
interval1 <- -qnorm((1-0.9)/2) # 90% multiplier
interval2 <- -qnorm((1-0.95)/2) # 95% multiplier
# Plot
zp1 <- ggplot(allModelFrame, aes(colour = modelName))
zp1 <- zp1 + geom_hline(yintercept = 0, colour = gray(1/2), lty = 2)
zp1 <- zp1 + geom_linerange(aes(x = Variable, ymin = Coefficient - SE*interval1,
ymax = Coefficient + SE*interval1),
lwd = 1, position = position_dodge(width = 1/2))
zp1 <- zp1 + geom_pointrange(aes(x = Variable, y = Coefficient, ymin = Coefficient - SE*interval2,
ymax = Coefficient + SE*interval2),
lwd = 1/2, position = position_dodge(width = 1/2),
shape = 21, fill = "WHITE")
zp1 <- zp1 + coord_flip() + theme_bw()
zp1 <- zp1 + ggtitle("Comparing several models")
print(zp1) # The trick to these is position_dodge().
@dsparks
Owner
dsparks commented Dec 18, 2012

Resulting image

@vinayakh
vinayakh commented Oct 2, 2013

The link for the data is not working. Can you please make the data available on some site to try out the scripts ? you can mail it at at gmail The blog is.R is very informative.

@bfoste01

This is fantastic! I needed something to make coefficient plots for lme4 models where I utilized imputation methods and analyzed over imputed data sets in parallel. This does the trick. Thank you.

@llthurman

This is fantastic, thank you. The problem I was previously having was that coefplot {arm} cannot handle lme or nlme objects and I seem to have some issues installing the coefplot2 package in my most recent RStudio version, so I went back to trusty ggplot2. I thought I would share some code that I modified based on your example to make it work for lme....

Fit several [lme] models with the same DV

now including random intercept of "Subject"

model1 <- lme(pid7 ~ ideo7 + female + age + south, random=~1|Subject, data = ANES)
model2 <- lm(pid7 ~ ideo7 + female + age + female:age, random=~1|Subject, data = ANES)
model3 <- lm(pid7 ~ ideo7, random=~1|Subject, data = ANES) # These are just arbitrary examples

Put model estimates into temporary data.frames

The only difference being that you have to specify the table component of the lme summary (tTable instead of coef)

model1Frame <- data.frame(Variable = rownames(summary(model1)$tTable),
Coefficient = summary(model1)$tTable[, 1],
SE = summary(model1)$coef[, 2],
modelName = "South Indicator")
model2Frame <- data.frame(Variable = rownames(summary(model2)$tTable),
Coefficient = summary(model2)$tTable[, 1],
SE = summary(model2)$tTable[, 2],
modelName = "Age Interaction")
model3Frame <- data.frame(Variable = rownames(summary(model3)$tTable),
Coefficient = summary(model3)$tTable[, 1],
SE = summary(model3)$tTable[, 2],
modelName = "Univariate")

everything else is the same...

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.