These data support Harrison et al (2017) Best Practice for Mixed Effects Models in Ecology and Evolution. This document contains all the code necessary to reproduce the simulations and Figure 1

Random Intercept Model Data Generation

#Set Seed for Reproducibility
  set.seed(5072917)

#Global Intercept 
  mu<-50

#Sample Size
  nsample<-100
  ngroups<-5
  n<-nsample*ngroups

#Global Slope
  beta1<-1

#Slope Measurements 
  body<-rnorm(n,100,20)
  
#Draw Global Intercepts
  alphas<-rnorm(ngroups,mu,30)

#Indicator for Population
  groupid<-sample(rep(1:ngroups,each=nsample))

#Linear Predictor 
  y<-rnorm(n,alphas[groupid] + beta1 * body,10)

#Make Data Frame
  ri.data<-data.frame(y=y,body=body,group=factor(groupid))
  


####### Predictions Using Simulated Values 
  
  #Make a Covariate Vector to Predict for
  body.pred<-seq(0,max(body),len=100)
  
  #Assemble Predictions for Each Intercept Value
  newdata<-data.frame(pred=unlist(lapply(alphas,function(x){x + beta1 *body.pred })))
  newdata$group<-factor(rep(seq(5),each=100))  
  newdata$alpha<-rep(alphas,each=100)  
  newdata$body<-rep(body.pred,5)
  
  #Dataframe for Global Intercept and Slope
  globaldata<-data.frame(pred=mu + beta1 * body.pred)
  globaldata$body<-body.pred
  
  #Subset Dataframes - this code is necessary to differentiate between the range of observed data (solid line in plot) and tracking fitted lines abck to the intercept (dashed lines)
  newdata.intjoin<-subset(newdata,body<42)
  newdata<-subset(newdata,body>42)
  globaldata.intjoin<-subset(globaldata,body<42)
  globaldata<-subset(globaldata,body>42)

Random Slope Model Data Generation

#Set RNG Seed
  set.seed(5072016)
  
#Uncorrelated Random Slopes with Global Beta = 1
  betas<-rnorm(ngroups,beta1,0.8)
  
#ReDraw Response Values 
  y.slope<-rnorm(n,alphas[groupid] + betas[groupid] * body,20)
  
#Make Data Frame
  rs.data<-data.frame(y=y.slope,body=body,group=factor(groupid))
  
##### Plot Using Just Predicted values
  m2coefs<-data.frame(alpha=alphas,beta=betas)
  
#Assemble Predictions for Each Intercept Value
  library(plyr)
  slope.preds<-alply(as.matrix(m2coefs),1,function(x) x[1] + x[2]*body.pred)
  
  newdata.slope<-data.frame(pred=unlist(slope.preds))
  newdata.slope$group<-factor(rep(seq(5),each=100))  
  newdata.slope$alpha<-rep(alphas,each=100)  
  newdata.slope$body<-rep(body.pred,5)
  
#Dataframe for Global Intercept and Slope
  globalintercept.slope<-mu 
  globalslope.slope<-beta1
  
globaldata.slope<-data.frame(pred=globalintercept.slope + globalslope.slope * body.pred)
  globaldata.slope$body<-body.pred

#Subset Dataframes -this code is necessary to differentiate between the range of observed data (solid line in plot) and tracking fitted lines abck to the intercept (dashed lines)
  newdata.slope.intjoin<-subset(newdata.slope,body<42)
  newdata.slope<-subset(newdata.slope,body>42)
  globaldata.slope.intjoin<-subset(globaldata.slope,body<42)
  globaldata.slope<-subset(globaldata.slope,body>42)

Plotting

library(ggplot2)

  #Choose Some Group Colours
  library(RColorBrewer)
  groupcols<-sample(brewer.pal(5,"Set2"))
  
  text1<-paste("y[i] == alpha[j] + beta * x[i]")
  intercepts<-c("alpha[1]","alpha[2]","alpha[3]","alpha[4]","alpha[5]")
  alphas.adjust<-alphas
  alphas.adjust[1]<-alphas[1]-5


####### RANDOM INTERCEPT PLOT

 #Lines and Intercept Cartoons
  g1<-ggplot(globaldata,aes(x=body,y=pred)) + geom_line(size=2) + ylim(c(-20,350)) + xlim(c(-20,max(body.pred))) + ylab("Dependent Variable y") + xlab("Predictor Variable x")
  g1 <- g1 +  geom_point(data=ri.data,shape=21,colour="black",aes(x=body,y=y,bg=group),size=3,alpha=0.4) + scale_fill_manual(values=groupcols)
  g1<- g1 + geom_line(data=newdata,aes(x=body,y=pred,colour=group),size=1.5)  + theme_bw() + scale_colour_manual(labels=seq(5),values=c(groupcols))
  g1<- g1 + geom_line(data=newdata.intjoin,aes(x=body,y=pred,colour=group),size=1.5,linetype=2,alpha=0.5)
  g1 <- g1 + geom_line(data=globaldata.intjoin,aes(x=body,y=pred),size=2,linetype=2,colour="black",alpha=0.4)
  g1<- g1+  theme(axis.text=element_blank(),axis.ticks=element_blank(),axis.title=element_text(size=18),legend.position="none",panel.grid.major = element_blank(),panel.grid.minor = element_blank())
  g1 <-g1 + annotate("text",x=30,y=290,label="Random Intercepts",size=6)
  g1 <- g1+ annotate("text",x=30,y=260,label=text1,parse=T,size=8)
  g1<- g1 + annotate("text",x=rep(-5,5),y=alphas.adjust,label=intercepts,parse=T,size=5) + annotate("text",x=-5,y=mu-4,label=paste("mu[group]"),parse=T,size=6)
  
####### RANDOM SLOPE PLOT
  
  text2<-paste("y[i] == alpha[j] + beta[j] * x[i]")
  intercepts<-c("alpha[1]","alpha[2]","alpha[3]","alpha[4]","alpha[5]")
  alphas.adjust<-alphas
  alphas.adjust[1]<-alphas[1]-5
  
    #Lines and Intercept Cartoons
  g2<-ggplot(globaldata.slope,aes(x=body,y=pred)) + geom_line(size=2) + ylim(c(-20,350)) + xlim(c(-20,max(body.pred))) + ylab("Dependent Variable y") + xlab("Predictor Variable x")
  g2 <- g2 +  geom_point(data=rs.data,shape=21,colour="black",aes(x=body,y=y,bg=group),size=3,alpha=0.4) + scale_fill_manual(values=groupcols)
  g2<- g2 + geom_line(data=newdata.slope,aes(x=body,y=pred,colour=group),size=1.5)  + theme_bw() + scale_colour_manual(labels=seq(5),values=c(groupcols))
  g2<- g2 + geom_line(data=newdata.slope.intjoin,aes(x=body,y=pred,colour=group),size=1.5,linetype=2,alpha=0.5)
  g2 <- g2 + geom_line(data=globaldata.slope.intjoin,aes(x=body,y=pred),size=2,linetype=2,colour="black",alpha=0.4)
  g2<- g2+  theme(axis.text=element_blank(),axis.ticks=element_blank(),axis.title=element_text(size=15),legend.position="none",panel.grid.major = element_blank(),panel.grid.minor = element_blank())
  g2 <-g2 + annotate("text",x=30,y=290,label="Random Intercepts and Slopes",size=6)
  g2 <- g2+ annotate("text",x=25,y=260,label=text2,parse=T,size=8)
  g2 <- g2 + annotate("text",x=rep(-5,5),y=alphas.adjust,label=intercepts,parse=T,size=5) + annotate("text",x=-5,y=mu-4,label=paste("mu[group]"),parse=T,size=6)
 
  
####### COMBINED PLOT 
  
  library(cowplot)
## 
## Attaching package: 'cowplot'
## The following object is masked from 'package:ggplot2':
## 
##     ggsave
  #Plot
    g3<-plot_grid(g1,g2,labels=c("A","B"),label_size=25)
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 5 rows containing missing values (geom_path).
    g3

  #Store Plots
    #ggsave('Mixed Effects Types.pdf',g3,width=15,height=6)
    #ggsave('Mixed Effects Types.tiff',g3,width=15,height=6)