Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
#############################################################
## 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()
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.