ploting_stimates_in_base.R

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


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 = 100)
p1()
dev.off()

png("p2.png",width=5, height = 5, units="in", res = 100)
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"))