##Here we present some simple mixed model examples in a simulated dataset

First we will set a seed for reproducibility. We have selected this one because it produces a clear plot. Feel free to change the seed and see how the outcome changes

set.seed(20)

We now load the packages we need

#libraries
library(nlme)
library(lme4)

Now we simulate the data we will fit models too

#random effect groups
clones<-6
#number of samples per random effect group
perclone<-20
clone<-gl(clones,perclone)
#categorical explanatory variable
treat<-rep(gl(2,perclone/2),clones)
treat<-as.numeric(treat)-1
#continuous explanatory variable (covariate)
treat.cov<-runif(clones*perclone,0,1)

#2 response variables - one for the categorical explanatory variable and one for the covariate
y.mean<-rep(rnorm(clones,0,1),each=perclone)
y.slope<-rep(rnorm(clones,0,0.5),each=perclone)
y.obs<-rnorm(perclone*clones,y.mean+treat+y.slope*treat,0.5)
y.obs.cov<-y.obs-y.slope*treat-treat+(1+y.slope)*treat.cov

We now fit the random intercepts versions of the models using nlme

#random intercept, treatment model
m1.rint<-lme(y.obs~treat,random=~1|clone, control = lmeControl(opt = "optim"))
summary(m1.rint)
## Linear mixed-effects model fit by REML
##  Data: NULL 
##        AIC      BIC    logLik
##   224.9004 235.9831 -108.4502
## 
## Random effects:
##  Formula: ~1 | clone
##         (Intercept)  Residual
## StdDev:    1.172938 0.5316437
## 
## Fixed effects: y.obs ~ treat 
##                  Value Std.Error  DF   t-value p-value
## (Intercept) -0.2931516 0.4837438 113 -0.606006  0.5457
## treat        1.1407645 0.0970644 113 11.752654  0.0000
##  Correlation: 
##       (Intr)
## treat -0.1  
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -2.3794186 -0.6752965  0.1082650  0.6576496  2.3083782 
## 
## Number of Observations: 120
## Number of Groups: 6
#random intercept covariate model
m1.rintcov<-lme(y.obs.cov~treat.cov,random=~1|clone, control = lmeControl(opt = "optim"))
summary(m1.rintcov)
## Linear mixed-effects model fit by REML
##  Data: NULL 
##        AIC      BIC    logLik
##   208.2389 219.3217 -100.1195
## 
## Random effects:
##  Formula: ~1 | clone
##         (Intercept)  Residual
## StdDev:    1.174101 0.4961598
## 
## Fixed effects: y.obs.cov ~ treat.cov 
##                  Value Std.Error  DF   t-value p-value
## (Intercept) -0.0355387 0.4874938 113 -0.072901  0.9420
## treat.cov    0.6133837 0.1534115 113  3.998290  0.0001
##  Correlation: 
##           (Intr)
## treat.cov -0.157
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -2.33587731 -0.64669678  0.02352256  0.61664940  2.05224191 
## 
## Number of Observations: 120
## Number of Groups: 6

We now fit the random slopes versions of the models using nlme

#random slope, treatment model
m1.rslope<-lme(y.obs~treat,random=~treat|clone, control = lmeControl(opt = "optim"))
summary(m1.rslope)
## Linear mixed-effects model fit by REML
##  Data: NULL 
##        AIC      BIC    logLik
##   217.0846 233.7087 -102.5423
## 
## Random effects:
##  Formula: ~treat | clone
##  Structure: General positive-definite, Log-Cholesky parametrization
##             StdDev    Corr  
## (Intercept) 1.0720340 (Intr)
## treat       0.4404704 0.382 
## Residual    0.4896137       
## 
## Fixed effects: y.obs ~ treat 
##                  Value Std.Error  DF   t-value p-value
## (Intercept) -0.2931516 0.4421970 113 -0.662943  0.5087
## treat        1.1407645 0.2008144 113  5.680692  0.0000
##  Correlation: 
##       (Intr)
## treat 0.293 
## 
## Standardized Within-Group Residuals:
##          Min           Q1          Med           Q3          Max 
## -2.249516136 -0.673480190  0.004823643  0.703378061  2.114185365 
## 
## Number of Observations: 120
## Number of Groups: 6
ranef(m1.rslope)
##   (Intercept)      treat
## 1   0.4725388  0.4082560
## 2   0.2751279 -0.4549530
## 3  -0.1686910  0.3879957
## 4  -0.7709775 -0.4163447
## 5   1.6298310  0.2514858
## 6  -1.4378293 -0.1764398
fixef(m1.rslope)
## (Intercept)       treat 
##  -0.2931516   1.1407645
#random slope covariate model
m1.rslopecov<-lme(y.obs.cov~treat.cov,random=~treat.cov|clone, control = lmeControl(opt = "optim"))
summary(m1.rslopecov)
## Linear mixed-effects model fit by REML
##  Data: NULL 
##        AIC      BIC    logLik
##   209.6554 226.2795 -98.82768
## 
## Random effects:
##  Formula: ~treat.cov | clone
##  Structure: General positive-definite, Log-Cholesky parametrization
##             StdDev    Corr  
## (Intercept) 1.1079246 (Intr)
## treat.cov   0.4116194 0.249 
## Residual    0.4827630       
## 
## Fixed effects: y.obs.cov ~ treat.cov 
##                  Value Std.Error  DF    t-value p-value
## (Intercept) -0.0369545 0.4606034 113 -0.0802305  0.9362
## treat.cov    0.6107964 0.2251688 113  2.7126153  0.0077
##  Correlation: 
##           (Intr)
## treat.cov 0.075 
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -2.16873970 -0.66078372  0.05075708  0.75091755  2.09459239 
## 
## Number of Observations: 120
## Number of Groups: 6
ranef(m1.rslopecov)
##   (Intercept)    treat.cov
## 1   0.5242008  0.295154222
## 2   0.2215720 -0.340176948
## 3  -0.1441830  0.254298446
## 4  -0.7619135 -0.431431787
## 5   1.6704060  0.219866382
## 6  -1.5100822  0.002289684
fixef(m1.rslopecov)
## (Intercept)   treat.cov 
## -0.03695446  0.61079645

We now plot the data and our model predictions to demonstrate the difference between the random intercepts and random slopes models

#set colour scheme for plots
cols<-c("#332288", "#88CCEE", "#117733", "#DDCC77", "#CC6677","#AA4499")

#plotting
par(mfrow=c(2,2),mar=c(3,3,1,1))

plot(y.obs~treat,type="n",xlab="",ylab="",xlim=c(-0.05,1.05),xaxt="n",las=1)
axis(side=1,at=c(0,1))
mtext(side=1,line=2,"Treatment")
mtext(side=2,line=2,"Response")
for(i in 1:clones){
  lines(c(0,1),c(fixef(m1.rint)[1]+ranef(m1.rint)[i,1],fixef(m1.rint)[1]+ranef(m1.rint)[i,1]+fixef(m1.rint)[2]*1),lwd=2,col=cols[i])}
for(i in 1:clones){
  treat.plot<-treat+i/50-0.06
  points(y.obs~treat.plot,subset=(clone==i),pch=21,cex=1.5,bg=cols[i])}
lines(c(0,1),c(fixef(m1.rint)[1],fixef(m1.rint)[1]+fixef(m1.rint)[2]*1),lwd=3,col="black",lty=2)

plot(y.obs~treat,type="n",xlab="",ylab="",xlim=c(-0.05,1.05),xaxt="n",las=1)
axis(side=1,at=c(0,1))
mtext(side=1,line=2,"Treatment")
mtext(side=2,line=2,"Response")
for(i in 1:clones){
  lines(c(0,1),c(fixef(m1.rslope)[1]+ranef(m1.rslope)[i,1],fixef(m1.rslope)[1]+ranef(m1.rslope)[i,1]+(fixef(m1.rslope)[2]+ranef(m1.rslope)[i,2])*1),lwd=2,col=cols[i])}
for(i in 1:clones){
  treat.plot<-treat+i/50-0.06
  points(y.obs~treat.plot,subset=(clone==i),pch=21,cex=1.5,bg=cols[i])}
lines(c(0,1),c(fixef(m1.rslope)[1],fixef(m1.rslope)[1]+fixef(m1.rslope)[2]*1),lwd=3,col="black",lty=2)

plot(y.obs.cov~treat.cov,type="n",xlab="",ylab="",xlim=c(-0.05,1.05),las=1)
mtext(side=1,line=2,"Covariate")
mtext(side=2,line=2,"Response")
treat.plot<-treat.cov
for(i in 1:clones){
  lines(c(0,1),c(fixef(m1.rintcov)[1]+ranef(m1.rintcov)[i,1],fixef(m1.rintcov)[1]+ranef(m1.rintcov)[i,1]+fixef(m1.rintcov)[2]*1),lwd=2,col=cols[i])}
for(i in 1:clones){
  points(y.obs.cov~treat.plot,subset=(clone==i),pch=21,cex=1.5,bg=cols[i])}
lines(c(0,1),c(fixef(m1.rintcov)[1],fixef(m1.rintcov)[1]+fixef(m1.rintcov)[2]*1),lwd=3,col="black",lty=2)

plot(y.obs.cov~treat.cov,type="n",xlab="",ylab="",xlim=c(-0.05,1.05),las=1)
mtext(side=1,line=2,"Covariate")
mtext(side=2,line=2,"Response")
for(i in 1:clones){
  lines(c(0,1),c(fixef(m1.rslopecov)[1]+ranef(m1.rslopecov)[i,1],fixef(m1.rslopecov)[1]+ranef(m1.rslopecov)[i,1]+(fixef(m1.rslopecov)[2]+ranef(m1.rslopecov)[i,2])*1),lwd=2,col=cols[i])}
for(i in 1:clones){
  treat.plot<-treat.cov+i/50-0.06
  points(y.obs.cov~treat.plot,subset=(clone==i),pch=21,cex=1.5,bg=cols[i])}
lines(c(0,1),c(fixef(m1.rslopecov)[1],fixef(m1.rslopecov)[1]+fixef(m1.rslopecov)[2]*1),lwd=3,col="black",lty=2)