Plotting model estimates in base R graphics

There are many ways of displaying data and graphs in R, here I propose a base graph solution to plot estimates and confidence intervals using a Generalized Linear Mixed Effect Model, but this code could be used for other kind of models. So, we have a Generalized Linear Effect Model and we want to plot the model estimates and its confidence interval. First we install needed packages and generate some data. Imagine we have a between subjects experiment with three conditions in which participants have to answer 30 trials as correct or incorrect and we want to know if the condition changes the number of correct answers.

This post takes this  webpage  as base. Check it out!

FULL CODE GIST 

if (!require("pacman")) install.packages("pacman")
pacman::p_load("lme4", "effects", "dplyr")

Generate data

parts = 30 # Our participants
trials = 30 # Our Trials
conds = 3 # How many conditions we have
response1 <- c();response2 <- c(); response3 <-c()
for(i in 1:parts){ 
  # As I'm not focusing on simulatng realistic data, the responses are totally trivial  
  response1 <- c(response1, rbinom(trials,1, 0.3+runif(1,-0.15, + 0.15) ))
  response2 <-c(response2, rbinom(trials,1, 0.4+runif(1,-0.25, + 0.25) ))
  response3 <-c(response3, rbinom(trials,1, 0.7+runif(1,-0.15, + 0.05) ))
  }
response <- c(response1, response2, response3)

condition <- c(rep(1, parts*trials), rep(2, parts*trials),rep(3, parts*trials))
participant <- sort(rep(1:(conds*parts), trials) )
trial_id <-rep(1:30, parts)
# our data frame, ready
df <- data.frame(participant, condition, response, trial_id)
df$condition <- as.factor(df$condition)

Modelling

Next step: fit a model with lme4.

m <- glmer(response ~  condition  + 
                  (1|participant) + (1|trial_id),
                family = binomial(), data =  df,
                na.action = na.omit, control=glmerControl(optimizer="bobyqa")) 

Create a fit and boundaries to plot

Now we have estimated the model. Is time to get the confindence intervals. There are several ways to get them. I'm using effects package because is quick and in this post I'm focusing on plotting rather than making models.

forplot.binom <- function(model, term){
  # This function takes the model and the term and gives the probabilites 
  # via the plogis() function we transform the effects output. Is better to have this in other script. 
  # you can modify it for getting estimates of non-binomial models, also.
  require(effects)
  forplot<-effect(term, model)$x
  forplot$fit<-plogis(effect(term, model)$fit)
  forplot$lower<-plogis(effect(term, model)$lower)
  forplot$upper<-plogis(effect(term, model)$upper)
  return(forplot)
}

(m.plot <- forplot.binom(m, "condition"))

Plot

Now we have the object m.plot which we will use for plotting. As I said before there are other ways to get the fit and bounds. You can use a different way and plot with the same script. You just have to be sure that the data.frame looks like the m.plot.

# Convert the m.plot into a thing to send to the plot 
fits = as.vector( m.plot$fit )
low = as.vector( m.plot$upper )
up = as.vector( m.plot$lower )
# Another helper function done by Eric-Jan Wagenmakers and Quentin F. Gronau
# in which all the plot is based 
# and taken from http://shinyapps.org/apps/RGraphCompendium/index.php
plotsegraph <- function(loc, fit, up, low,  color = "grey", linewidth = 2) {
segments(x0 = loc, x1 = loc, y0 = fit  , y1 = up , col = color, 
          lwd = linewidth)
segments(x0 = loc, x1 = loc, y0 = fit  , y1 = low, col = color,
          lwd = linewidth)
}
# A nice plotting way is make the plot as a function. Is a easy way to save the 
# plot in different formats (look at the end ;).

p1 <- function(){
par(cex.main = 1.5, mar = c(5, 6, 4, 5) + 0.1, mgp = c(3.5, 1, 0), cex.lab = 1.5, 
    font.lab = 2, cex.axis = 1.3, bty = "n", las = 1)
x <- c(1, 2, 3) # How many conditions we have
# The logic following this plot is creating first an empty canvas and plotting y outside of the limits defined by
# ylim and xlim. 
plot(x, y = c(-10, -10, -10), type = "p", ylab = "", xlab = " ", cex = 1.5, 
     ylim = c(0.1, 0.9), xlim = c(1, 3.5), lwd = 2, pch = 5, axes = F, main = "Title: ...")
# Then we make the axis labels. 
axis(1, at = c(1.5, 2.5, 3.5), labels = c("One group", "Other", "Another one"))
seq_tick = c(seq(0,1, 0.1) ) # change values in seq for different labels 
axis(2, pos = 1,labels = seq_tick, at = seq_tick  ) 
mtext("Here come the conditions", side = 1, line = 3, cex = 1.5, font = 2)
 
par(las = 0)
mtext("Probability of correct item", side = 2, line = 2, cex = 1.5, font = 2)
x <- c(1.5, 2.5, 3.5)
# Plotting the fits
points(x, c(fits), cex = 1.5, lwd = 5, pch = 19)
# Plotting the boundaries with the helper function
plot.errbars <- plotsegraph(loc=x, 
                            fit=fits, 
                            low=low, 
                            up=up,
                            4, color = "black") 
# I wanted to add a line marking the 0.5 
lines(c(1.2, 3.7), c(0.5,0.5) , lty = 2,lwd=2)
#  a grid : )
for (i in c(seq(0.1, 0.9,.2))){
   
  lines(c(1.2,3.7), c(i,i), lwd=0.3, lty=2)
}
# and some text
text(x= 3, y = 0.51,label = "A text", adj=c(0,0))

}
# execute the next line to see the plot. 
p1()

Creating a simple plotting function

It looks good, but what about making the function completely reusable?

plot_estimates <- function(m.plot, linewidth=4, grid = TRUE, 
                           ylims=c(0,1), plot_line_at=NULL, label_conds = NULL, title_plot = NULL, y_label = NULL){

fits = as.vector( m.plot[,2] )
low = as.vector( m.plot[,3] )
up = as.vector( m.plot[,4] )

if(is.null(label_conds)) {label_conds = as.vector( m.plot[,1])}
if(is.null(title_plot)) {title_plot = ""}
if(is.null(y_label)) {y_label = ""}

plotsegraph <- function(loc, fit, up, low,  color = "grey", linewidth = 2) {
segments(x0 = loc, x1 = loc, y0 = fit  , y1 = up , col = color, 
          lwd = linewidth)
segments(x0 = loc, x1 = loc, y0 = fit  , y1 = low, col = color,
          lwd = linewidth)
}
number_of_conds = length(fits)
par(cex.main = 1.5, mar = c(5, 6, 4, 5) + 0.1, mgp = c(3.5, 1, 0), cex.lab = 1.5, 
    font.lab = 2, cex.axis = 1.3, bty = "n", las = 1)
x1 <- c(1:number_of_conds) # How many conditions we have
plot(x1, c(rep(-10, number_of_conds)), type = "p", ylab = "", xlab = "", cex = 1.5, 
     ylim = c(0,1), xlim = c(1, number_of_conds+0.5), lwd = 2, pch = 5, axes = F, main = title_plot)

axis(1, at = c(1:number_of_conds)+0.5 , labels = label_conds)
axis(2, pos = 1, labels = c(seq(0,1, 0.1) ), at = seq(0,1, 0.1)  )
mtext("Here come the conditions", side = 1, line = 3, cex = 1.5, font = 2)
 
par(las = 0)
mtext(y_label, side = 2, line = 2, cex = 1.5, font = 2)
x2 <- c(1:number_of_conds) + .5
points(x2, c(fits), cex = 1.5, lwd = linewidth, pch = 19)
plot.errbars <- plotsegraph(loc=x2, 
                            fit=fits, 
                            low=low, 
                            up=up,
                            linewidth, color = "black") 
if(!is.null(plot_line_at)) lines(c(1, number_of_conds+.5), c(plot_line_at,plot_line_at) , lty = 2,lwd=2)

if(grid){
  for (i in c(seq(0.1, 0.9,.2))){
    lines(c(1,number_of_conds+.5), c(i,i), lwd=0.3, lty=2)
    }}

}

# we fit a model with less conditions
m2 <-  glmer(response ~  condition  + 
                  (1|participant) + (1|trial_id),
                family = binomial(), data =  df[which(df$condition != 3),],
                na.action = na.omit, control=glmerControl(optimizer="bobyqa")) 
m2.plot <- forplot.binom(m2, "condition")

plot_estimates(m2.plot, grid = TRUE, ylims=c(0.1,0.9), plot_line_at=0.5, label_conds = c("One group", "Another"), title_plot = "Which is the best group?", y_label = c("Mean estimated score"))

Saving plots

# Then if you like it, you can save it
tiff("p1.tiff", width=5, height = 5, units="in", res = 100)
p1()
dev.off()
# But rather than use the tiff, you would prefer make a png to put on a slide: 
png("p1.png", width=5, height = 5, units="in", res = 160)
p1()
dev.off()