Create a gist now

Instantly share code, notes, and snippets.

What would you like to do?
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().
Owner

dsparks commented Dec 18, 2012

Resulting image

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.

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 commented Nov 4, 2015 edited

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)$tTable[, 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...

Gavago commented Jun 13, 2017 edited

A wonderful example. I added position_dodge() to my existing code to make the comparisons I wanted.

I have a question though regarding the order of the colour factor levels...

My data frame coefs:

     Estimate         se     lowerCI     upperCI  dyads            vars
3   0.106478735 0.06169716 -0.01444769  0.22740516    All         Kinship
2  -2.591911491 0.20060472 -2.98509674 -2.19872624    All  Age difference
4  -0.169305922 0.07252870 -0.31146217 -0.02714967    All Rank difference
31  0.007409558 0.10206631 -0.19264041  0.20745952 Female         Kinship
21 -1.944700523 0.25320891 -2.44099000 -1.44841105 Female  Age difference
41 -0.220988389 0.10169709 -0.42031468 -0.02166210 Female Rank difference
32  0.137868642 0.06934290  0.00195656  0.27378072   Male         Kinship
22 -2.880896180 0.24594286 -3.36294418 -2.39884818   Male  Age difference
42 -0.136386916 0.08218800 -0.29747540  0.02470157   Male Rank difference

My code:

g<-ggplot(coefs, aes(vars, Estimate, colour=dyads)) + 
  theme(axis.title.x =element_text(margin=margin(0,20,0,0)))+
  theme(axis.title.x = element_blank())+
  
  geom_hline(yintercept=0, lty=2, lwd=1, colour="grey50") +
  
  geom_errorbar(aes(ymin=lowerCI, ymax=upperCI), 
                lwd=1, width=0, position = position_dodge(width = 0.5)) + 
  labs(x ="")+
  geom_point(size=4, pch=21,position = position_dodge(width = 0.5)) +
  
  theme(panel.background = element_rect(fill = "white"))
g

When I add coord_flip() at the end of this code, the lines are no longer ordered by the levels of the factor coefs$dyads. As you can see in your example also, @dsparks,South Indicator is the first level of the factor modelName, but its lines are below the lines of Age Interaction and Univariate. Do you know how to custom rearrange these lines such that they appear in the same order as the levels of the colour factor?

I realize I can simply omit coord_flip(), save the figure, and rotate it by hand, but rotating the coordinates beforehand is, of course, much more desirable.

Thanks for any insight.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment