Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
###"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"))
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.