ggplot2 Caterpillar Plot

ggplot2_caterpillar.R

# Load libraries
lapply(c("R2MLwiN","arm","ggplot2","RColorBrewer","grid"),library, character.only=TRUE)

# Set working directory and load data
setwd("~/R")
data(tutorial)
str(tutorial)

# Fit a linear varying intercept multilevel model with group level predictors of normexam ~ sex + standlrt + vrband + schgend + avslrt + schav + (1|school)
fit.1 <-lmer(normexam ~ sex + standlrt + vrband + schgend + avslrt + schav + (1|school), data=tutorial)
display(fit.1)

# Display "Random"/"Varying" i.e. school effects
ranef(fit.1)

# Extract higher level residuals
ranef <-ranef(fit.1)
ranef.se<-se.ranef(fit.1)
u0mn <-ranef$school
u0se <-ranef.se$school
# Rank residuals
u0rankmn = rank(u0mn)
u0hi <- u0mn + (1.96*u0se)
u0lo <- u0mn  - (1.96*u0se)

# Combine into data.frame
d <-data.frame(u0mn,u0rankmn, u0hi, u0lo)

# Set colors
palette <- brewer.pal("Greys", n=9)
color.background = palette[1]
color.grid.major = palette[5]
color.axis.text = palette[6]
color.axis.title = palette[7]
color.title = palette[9]

# Plot
k <- ggplot()+
  geom_pointrange(data=d,mapping=aes(x=d$u0rankmn, y=d$X.Intercept., ymin=d$X.Intercept..2,ymax=d$X.Intercept..1), position="identity", width=0.5, size=1, color="black")+
  theme_bw() +
  # Set the entire chart region to a light gray color
  theme(panel.background=element_rect(fill=color.background, color=color.background)) +
  theme(plot.background=element_rect(fill=color.background, color=color.background)) +
  theme(panel.border=element_blank(), axis.line=element_line(), axis.line.y=element_blank()) +  
  # Format the grid
  theme(panel.grid.major=element_line(colour=color.background,size=.75)) +
  theme(axis.ticks=element_blank()) +
  theme(axis.line.x =element_line(colour="#535353", size=.75))+
  theme(axis.line.y =element_line(colour="#535353", size=.75))+
  # Dispose of the legend
  theme(legend.position="none") +
  # Set title and axis labels, and format these and tick marks
  ggtitle("Residuals: Intercept")+ 
  theme(plot.title=element_text(face="bold",hjust=-.08,vjust=2,colour="#3C3C3C",size=20)) +
  xlab("Rank of residuals") +
  ylab("Intercept") +
  theme(axis.text.x=element_text(size=17,colour="#535353",face="bold")) +
  theme(axis.text.y=element_text(size=17,colour="#535353",face="bold")) +
  theme(axis.title.y=element_text(size=17,colour="#535353",face="bold",vjust=1.5)) +
  theme(axis.title.x=element_text(size=17,colour="#535353",face="bold",vjust=-.5)) +
  # Big bold line at y=0
  geom_hline(yintercept=0,size=1.2, alpha=0.7,colour="#EF3B2C", linetype="twodash")+
  # Plot margins and finally line annotations
  theme(plot.margin = unit(c(1, 1, .5, .7), "cm"))