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