Loading libraries

If you haven’t already done it, you’ll need to install the ggplot2 library and all its dependencies before running the next command.

library(ggplot2)    
## Warning: package 'ggplot2' was built under R version 3.3.2
library(lme4)   
## Loading required package: Matrix
library(boot)    # will be used later for the 'inv.logit' function

First example

Run a model, generate predicted values and model uncertainty (error bars) and plot. A thorough understanding of ggplot is beyond what I can show here, but here are a couple of examples that may prove useful. The first is a plot of myr1 from the simple predictors example, so we’ll recreate that model here.

Recall that for this data set, there was only one delay ratio - the SS delay was 2.5 and the LL delay was 7.5; only the SS and LL magnitudes varied.

myd<-read.csv("Example simple predictors VG.dat", na.string="")  
myd<-subset(myd, !is.na(Waited))  
myd$logMagRatio <- log(myd$LLmagn/myd$SSmagn)
myd$logDelayRatio <- log((myd$LLdelay+1)/(myd$SSdelay+1)) 
myr1<-glmer(Waited~logMagRatio+logDelayRatio-1+(logMagRatio+logDelayRatio-1|Subject), data=myd, family=binomial)
## Warning: Some predictor variables are on very different scales: consider
## rescaling

Note, EVERY variable in the model being plotted must be included in the newdat grid as well as referenced in the ggplot commands below.

newdat <- expand.grid(
  logMagRatio = c(0, .5, 1, 1.5, 2, 2.5) 
  , logDelayRatio =c(.8873)
  , Waited = 0
)

modelToPlot<-myr1 

# The following code was obtained on-line to extract the correct error bars (plo/phi) and put them in newdat.
mm <- model.matrix(terms(modelToPlot),newdat)
newdat$Waited <- mm %*% fixef(modelToPlot)
pvar1 <- diag(mm %*% tcrossprod(vcov(modelToPlot),mm))
newdat <- data.frame(
  newdat
  , ProbWaited = inv.logit(newdat$Waited)
  , plo = inv.logit(newdat$Waited-sqrt(pvar1))
  , phi = inv.logit(newdat$Waited+sqrt(pvar1))
)

newdat   # if you want to look at what we're going to plot. 
##   logMagRatio logDelayRatio      Waited ProbWaited       plo       phi
## 1         0.0        0.8873 -1.15522403  0.2395362 0.2043098 0.2787088
## 2         0.5        0.8873 -0.78143398  0.3140109 0.2802839 0.3498237
## 3         1.0        0.8873 -0.40764392  0.3994772 0.3674329 0.4324058
## 4         1.5        0.8873 -0.03385386  0.4915373 0.4575810 0.5255719
## 5         2.0        0.8873  0.33993620  0.5841750 0.5440696 0.6231959
## 6         2.5        0.8873  0.71372626  0.6712240 0.6241725 0.7150739

Plotting behavior as a function of logMagRatio. Note, the LLdelay is three times longer than the SSdelay so it isn’t surprising that when the LLmagn = SSmagn (log(MagRatio) = 0) that people were only waiting 25% of the time

g1 <- ggplot(newdat, aes(x=logMagRatio, y= ProbWaited))+geom_line()+geom_point()+
    ylim(0,1)   
g1 + geom_errorbar(aes(ymin = plo, ymax = phi), linetype=1, width=.02)+
  xlab("Log(LLMagnitude/SSMagnitude)")+ylab("P(Waited)\n")+theme_bw()+theme(text = element_text(size=18))

New more complex model

Now, we’re going to run and plot myr2.1 from the simple predictions example.

myd$c.Sex<-ifelse(myd$Sex=="f", 1, -1)
myr2.1<-glmer(Waited~logMagRatio+logDelayRatio
   + c.Sex:logMagRatio+c.Sex:logDelayRatio-1+(logMagRatio+logDelayRatio-1|Subject), data=myd, family=binomial)
## Warning: Some predictor variables are on very different scales: consider
## rescaling
newdat <- expand.grid(
  logMagRatio = c(0, .5, 1, 1.5, 2, 2.5) 
  , logDelayRatio =c(.8873)
  , c.Sex = c(-1, 1)
  , Waited = 0
)

modelToPlot<-myr2.1 

mm <- model.matrix(terms(modelToPlot),newdat)
newdat$Waited <- mm %*% fixef(modelToPlot)
pvar1 <- diag(mm %*% tcrossprod(vcov(modelToPlot),mm))
newdat <- data.frame(
  newdat
  , ProbWaited = inv.logit(newdat$Waited)
  , plo = inv.logit(newdat$Waited-sqrt(pvar1))
  , phi = inv.logit(newdat$Waited+sqrt(pvar1))
)

newdat$Sex <- ifelse(newdat$c.Sex==1, "Female", "Male")
newdat$logDelay <- as.factor(newdat$logDelay)

# Plotting behavior as a function of MagRatio
g1 <- ggplot(newdat, aes(x=logMagRatio, y= ProbWaited, group=Sex, linetype=Sex))+geom_line()+geom_point()+
  ylim(0,1)  
g1 + geom_errorbar(aes(ymin = plo, ymax = phi), linetype=1, width=.02)+
  xlab("Log(LLMagnitude/SSMagnitude)")+ylab("P(Waited)\n")+theme_bw()+theme(text = element_text(size=18))

Repeat of above but using position_dodge to avoid error bar overlap.

pd <- position_dodge(width = 0.06)
g1 <- ggplot(newdat, aes(x=logMagRatio, y= ProbWaited, group=Sex, linetype=Sex))+geom_line(position = pd)+geom_point(position = pd)+
  ylim(0,1)
g1 + geom_errorbar(aes(ymin = plo, ymax = phi), linetype=1, width=.02, position = pd)+
  xlab("Log(LLMagnitude/SSMagnitude)")+ylab("P(Waited)\n")+theme_bw()+theme(text = element_text(size=18))

New plot involving single group

Here I’m running and plotting the results from myr1 from the Example one group.R file. Problem arises in creating newdat with the complex predictor names, so redoing the analysis.

myd<-read.csv("Example one group data glmer.dat")  
myd$logMagRatio <- log(myd$LLmagn/myd$SSmagn)
myd$logDelay <- log(myd$LLdelay+1)
myr1<-glmer(Waited~logMagRatio+logDelay-1+(logMagRatio+logDelay-1|Subject), data=myd, family=binomial)

newdat <- expand.grid(
  logMagRatio = seq(0,1,.025) 
  , logDelay =c(2, 3.5, 5)
  , Waited = 0
)

modelToPlot<-myr1 

mm <- model.matrix(terms(modelToPlot),newdat)
newdat$Waited <- mm %*% fixef(modelToPlot)
pvar1 <- diag(mm %*% tcrossprod(vcov(modelToPlot),mm))
newdat <- data.frame(
  newdat
  , ProbWaited = inv.logit(newdat$Waited)
  , plo = inv.logit(newdat$Waited-sqrt(pvar1))
  , phi = inv.logit(newdat$Waited+sqrt(pvar1))
)

# Plotting behavior as a function of logMagRatio.  Note that the beta for logMagRatio is much greater than
# that for logDelay.
pd <- position_dodge(width = 0.06)
g1 <- ggplot(newdat, aes(x=logMagRatio, y= ProbWaited, group=logDelay, colour=logDelay))+geom_line(position = pd)+geom_point(position = pd)+
  ylim(0,1)   
g1 + geom_errorbar(aes(ymin = plo, ymax = phi), linetype=1, width=.02, position = pd)+
  xlab("Log(LLMagnitude/SSMagnitude)")+ylab("P(Waited)\n")+theme_bw()+theme(text = element_text(size=18))
## Warning: position_dodge requires non-overlapping x intervals

## Warning: position_dodge requires non-overlapping x intervals

That is pretty messy, so let’s try plotting the same results in a couple of different ways.

Ribbon plot

I’m going to use error ribbons to show the variability so there’s no reason to dodge. In the ribbon syntax, the ‘fill’ is used for coloring the ribbons to match the lines, and alpha to create transparency. I eventually want to break out the plots using facet_grid, so I need to change logDelay into a factor.

newdat$logDelay <- factor(newdat$logDelay)
g2 <- ggplot(newdat, aes(x=logMagRatio, y= ProbWaited, group=logDelay)) + geom_line(aes(color=logDelay))+ ylim(0,1)   
g2 + geom_ribbon(aes(ymin = plo, ymax = phi, fill=logDelay), alpha=.3)+
  xlab("Log(LLMagnitude/SSMagnitude)")+ylab("P(Waited)\n")+theme_bw()+theme(text = element_text(size=18))

Back to error bars but plotted separately

g2 <- ggplot(newdat, aes(x=logMagRatio, y= ProbWaited))+geom_line()+geom_point()+
  ylim(0,1)   
g2 + geom_errorbar(aes(ymin = plo, ymax = phi), linetype=1, width=.02)+
  xlab("Log(LLMagnitude/SSMagnitude)")+ylab("P(Waited)\n")+theme_bw()+theme(text = element_text(size=18))+facet_grid(logDelay~.)

Same thing but using an error ribbon.

g2 <- ggplot(newdat, aes(x=logMagRatio, y= ProbWaited, group=logDelay))+geom_line(aes(color=logDelay))+ ylim(0,1)   
g2 + geom_ribbon(aes(ymin = plo, ymax = phi, fill=logDelay), alpha=.3)+
  xlab("Log(LLMagnitude/SSMagnitude)")+ylab("P(Waited)\n")+theme_bw()+theme(text = element_text(size=18))+facet_grid(logDelay~.)

Back to error bars but different x-axis and delay values along right axis.

I don’t think in terms of logs, so let’s create new columns with the backtransformed values.

newdat$MagRatio <- exp(newdat$logMagRatio)
newdat$Delay <- exp(as.numeric(as.character(newdat$logDelay))) - 1  # complicated because I changed it to a factor earlier
newdat$Delay <- as.factor(round(newdat$Delay))                      # removes ugly decimal places

g3 <- ggplot(newdat, aes(x=MagRatio, y= ProbWaited))+geom_line()+geom_point()+
  ylim(0,1)   
g3 + geom_errorbar(aes(ymin = plo, ymax = phi), linetype=1, width=.02)+
  xlab("LLMagnitude/SSMagnitude")+ylab("P(Waited)\n")+theme_bw()+theme(text = element_text(size=18))+facet_grid(Delay~.)

  • Updated: 9/22/23