ploting_stimates_in_R.R

#############################################################
## Estimating and plotting Generalized Linear Mixed Models ##
#############################################################

# This script will
# 1) Simulate some experimental data
# 2) Run a glmer with afex
# 3) Estimate predicted marginal means with emmeans
# 4) Plot estimates

# Install and load packages
if (!require("pacman")) install.packages("pacman")
pacman::p_load("ggplot2", "lme4", "emmeans","afex")

# Simulate some date by trial
parts = 30 # Our participants
trials = 60 # 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.35, + 0.35) ))
  response2 <-c(response2, rbinom(trials,1, 0.4+runif(1,-0.75, + 0.75) ))
  response3 <-c(response3, rbinom(trials,1, 0.7+runif(1,-0.15, + 0.05) ))
}
response <- c(response1, response2, response3)

condition <- c(rep("a", parts*trials), rep("b", parts*trials),rep("c", 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)
contrasts(df$condition) <- contr.sum

# Modelling with afex 
m <- mixed(response ~  condition  + 
             (1|participant) + (1|trial_id), method = "LRT",
           family = binomial(), data =  df) 
# Calculate estimated marginal means with emmeans. 
# In earlier scripts I used lsmeans but if you are using 
# lsmeans take into account that you have to fit the model
# with lme4 instead of afex 

means <- emmeans::emmeans(m, "condition", type = "response")
means 

# Helper function to extract the fits and confidence intervals
forplot.means <- function(m, effects){
  # @m = linear mixed effect model
  # @effects =  condition which we want to extract marginal means
  m.plot <- emmeans::emmeans(m, effects, type = "response", adjust = "holm")
  m.plot <- summary(m.plot)
  m.plot$lower <- m.plot$asymp.LCL
  m.plot$upper <- m.plot$asymp.UCL
  m.plot$fit <- m.plot$prob
  return(as.data.frame(m.plot) )
}
m.plot <- forplot.means(m, effects =  "condition")

# First we do a nice and clean plot in ggplot2 

theme_set(theme_bw())
my_theme <- theme(panel.grid.major.x = element_blank(),
                  panel.grid.minor.x= element_blank(),
                  panel.background = element_blank(),
                  panel.border = element_blank(),
                  strip.text.x = element_text(size=14), strip.text.y = element_text(size=16),
                  strip.background = element_rect(colour="black"),
                  axis.text.x  = element_text(color="black", size=14), 
                  axis.text.y  = element_text(color="black", size=14),
                  axis.title.x = element_text(face="bold", size=18),
                  axis.title.y = element_text(face="bold",size=18),
                  axis.ticks =  element_blank(),
                  #axis.ticks.y = element_blank(),
                  legend.title = element_text(size=14), legend.text = element_text(size = 14),
                  plot.title = element_text(vjust = 1.3, face="bold", size=20),
                  plot.margin = unit(c(1.5, 1.5, 1.5, 1.5), "cm") 
)

###############################
## ggplot2 version
###############################

pd <- position_dodge(0.1) # move them .05 to the left and right
plot.1 <- ggplot(m.plot, aes(x = condition, y = fit)) + 
  geom_linerange(aes(ymin=lower, ymax=upper), colour="black",  size = 2)+
  geom_point(position=pd, size=6) + 
  my_theme +
  scale_y_continuous(limits = c(0, 1), breaks = c(0, .25, .5, .75, 1)) + 
  ylab("Probability") + xlab("Condition") +
  scale_x_discrete(labels = c("Group A", "Group B", "Group C" ))

plot.1

ggsave("estimates_plot.png", plot = plot.1,   width = 7, height = 5)


###############################
## base version
###############################

fits <- c(as.vector(m.plot$fit)) 
up <- c(as.vector(m.plot$upper))
low <- c(as.vector(m.plot$lower) )

# helper function
plotsegraph <- function(loc, fit, up, low,  color = "grey", linewidth = 2) {
  #w <- wiskwidth/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)
}


p1 <- function(){
  par(cex.main = 1.5, mar = c(6, 5, 4, 4) + 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, 4)
  plot(x,rep(-10, length(x)), type = "p", ylab = "", xlab = " ", cex = 1.5, 
       ylim = c(0.15, 0.85), xlim = c(1, 4), lwd = 2, pch = 5, axes = F, main = "")
  # X axis labels 
  axis(1, at = c(seq(1.5, 3.5, 1)), labels = c("Group A", "Group B", "Group C"),tick = FALSE, lwd = 0)
  # Y axis labels
  axis(2, pos = 1,labels = gsub(as.character( seq(0,1, 0.1)), pattern = "0", replacement = ""), at = seq(0,1, 0.1)  )
  # X axis text
  mtext("Condition", side = 1, line = 4, cex = 1.5)
  
  par(las = 0)
  # Y axis
  mtext(expression(paste("Probability")), side = 2, line = 2, cex = 1.5, font = 2)
  x <- c(1.5, 2.5, 3.5)
  points(x, c(fits), cex = 1.5, lwd = 5, pch = 19)
  plot.errbars <- plotsegraph(loc=x, fit=fits, low=low, up=up,
                              4, color = "black") 
  lines(c(1.2, 6.7), c(0.5,0.5) , lty = 2,lwd=2)
  # Adding the grid
  for (i in c(seq(0.1, 0.9,.2))){
    
    lines(c(1.2,6.7), c(i,i), lwd=0.3, lty=2)
  }
  # add a chance level line
  text(x= 1.6, y = 0.45,label = "Chance level", adj=c(0,0))
  
}

# Plot the graphs 
p1()

# 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()