Create a gist now

Instantly share code, notes, and snippets.

Embed
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

This comment has been minimized.

Show comment
Hide comment
Owner

dsparks commented Dec 18, 2012

Resulting image

@vinayakh

This comment has been minimized.

Show comment
Hide comment
@vinayakh

vinayakh 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.

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 comment has been minimized.

Show comment
Hide comment
@bfoste01

bfoste01 Jul 10, 2015

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.

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 comment has been minimized.

Show comment
Hide comment
@llthurman

llthurman Nov 4, 2015

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...

llthurman commented Nov 4, 2015

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

This comment has been minimized.

Show comment
Hide comment
@Gavago

Gavago Jun 13, 2017

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.

Gavago commented Jun 13, 2017

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.

@superdayv

This comment has been minimized.

Show comment
Hide comment
@superdayv

superdayv Oct 28, 2017

Is it possible to do this with models that were fit using "brms"?

superdayv commented Oct 28, 2017

Is it possible to do this with models that were fit using "brms"?

@wraymond

This comment has been minimized.

Show comment
Hide comment
@wraymond

wraymond Jan 31, 2018

This is awesome! Thanks.

This is awesome! Thanks.

@stapial

This comment has been minimized.

Show comment
Hide comment
@stapial

stapial Mar 7, 2018

This is awesome, thanks!

Quick question, is there anyway to do the same but have each model in a different column and not stacked? I have 6 binary regression models (they don't have the same dependent variable). It looks messy when I lump them all together. I would prefer to have a column for each model with its respective coefficients in the same plot. Any ideas on how to do this???

Thanks!

stapial commented Mar 7, 2018

This is awesome, thanks!

Quick question, is there anyway to do the same but have each model in a different column and not stacked? I have 6 binary regression models (they don't have the same dependent variable). It looks messy when I lump them all together. I would prefer to have a column for each model with its respective coefficients in the same plot. Any ideas on how to do this???

Thanks!

@ksorby

This comment has been minimized.

Show comment
Hide comment
@ksorby

ksorby Mar 12, 2018

@stapial
You can use ggplot2's facet functions to split the models into columns:

zp1<- zp1 + facet_grid(.~modelName)

ksorby commented Mar 12, 2018

@stapial
You can use ggplot2's facet functions to split the models into columns:

zp1<- zp1 + facet_grid(.~modelName)

@stapial

This comment has been minimized.

Show comment
Hide comment
@stapial

stapial Mar 16, 2018

perfect! Thank you so much @ksorby!

stapial commented Mar 16, 2018

perfect! Thank you so much @ksorby!

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