Peril 4b: Problems caused by having exceptional slopes

Here we simulate an example of situations where there is only a relationship in one random effect group and demonstrate the power of random slopes models to deal with this

Load required R packages

library(lme4)
## Loading required package: Matrix

Here we define the parameters used in generating simulated datasets

#Number and size of groups used as random effects
groups<-10
groupsize<-20

#Calculate total sample size
n<-groups*groupsize

#Simulate values of an explanatory variable (independent from random effects)
bodysize<-rnorm(n,mean=0,sd=20)

#Simulate group IDs
group_ids<-rep(seq(1,groups,1),each=groupsize)


#Define parameters responsible for controlling values of the response variable
#residual variation
resid<-3
#Slopes in normal groups
slopeN<-0
#Define two potential slope values for the exception group
slope1<-0.1
slope2<-0.2

#Define standard deviation for the random effect hyperparameter (intergroup residual)
int.resid<-3

#Create a vector defining whether an observation is in the exception group or not
slopeN_v<-rep(0,n)
slope1_v<-c(rep(0,n-groupsize),rep(slope1,groupsize))
slope2_v<-c(rep(0,n-groupsize),rep(slope2,groupsize))

#Create empty vectors to store values of the response variable
signalsizeN<-rep(NA,n)
signalsize1<-rep(NA,n)
signalsize2<-rep(NA,n)

Here we create empty matrices to store results from the simulations

IoutputN<-matrix(NA,nr=100,nc=2)
Ioutput1<-matrix(NA,nr=100,nc=2)
Ioutput2<-matrix(NA,nr=100,nc=2)

SoutputN<-matrix(NA,nr=100,nc=2)
Soutput1<-matrix(NA,nr=100,nc=2)
Soutput2<-matrix(NA,nr=100,nc=2)

Now we can conduct the simulations

NB. We are suppressing the warning messages from model fitting from being printed here for presentation purposes. We would recommend re-running after deleting warning=FALSE and message=FALSE from the line below this to see what warnings arise and when they do if you are interested.

for(i in 1:100){
  
    #Simulate values of the response variable
    int.resids<-rep(rnorm(groups,0,int.resid),each=10)
    signalsizeN<-rnorm(n,slopeN_v*bodysize,resid)
    signalsize1<-rnorm(n,slope1_v*bodysize,resid)
    signalsize2<-rnorm(n,slope2_v*bodysize,resid)
    
    #Fit random intercepts models to:
    #modN - no exception slope
    #mod1 - exception slope with smaller effect size
    #mod2 - exception slope with larger effect size
    modN<-lmer(signalsizeN~bodysize+(1|group_ids))
    mod1<-lmer(signalsize1~bodysize+(1|group_ids))
    mod2<-lmer(signalsize2~bodysize+(1|group_ids))
    
    #Fit random slopes models to:
    #modN - no exception slope
    #mod1 - exception slope with smaller effect size
    #mod2 - exception slope with larger effect size
    modNa<-lmer(signalsizeN~bodysize+(bodysize|group_ids))
    mod1a<-lmer(signalsize1~bodysize+(bodysize|group_ids))
    mod2a<-lmer(signalsize2~bodysize+(bodysize|group_ids))
    
    #Store model estimates for both models in the results matrices
    IoutputN[i,1]<-fixef(modN)[2]
    Ioutput1[i,1]<-fixef(mod1)[2]
    Ioutput2[i,1]<-fixef(mod2)[2]
    IoutputN[i,2]<-summary(modN)[[10]][2,2]
    Ioutput1[i,2]<-summary(mod1)[[10]][2,2]
    Ioutput2[i,2]<-summary(mod2)[[10]][2,2]
    
    SoutputN[i,1]<-fixef(modNa)[2]
    Soutput1[i,1]<-fixef(mod1a)[2]
    Soutput2[i,1]<-fixef(mod2a)[2]
    SoutputN[i,2]<-summary(modNa)[[10]][2,2]
    Soutput1[i,2]<-summary(mod1a)[[10]][2,2]
    Soutput2[i,2]<-summary(mod2a)[[10]][2,2]

}

We now do some summarising a plotting of the results

First we work out the number of simulation runs (out of 100) in which we detected an overall positive slope for both the random intercepts model and random slopes model

#For random intercepts models
sum((IoutputN[,1]-IoutputN[,2]*1.96)>0) #null
## [1] 1
sum((Ioutput1[,1]-Ioutput1[,2]*1.96)>0) #smaller effect
## [1] 17
sum((Ioutput2[,1]-Ioutput2[,2]*1.96)>0) #larger effect
## [1] 32
#for random slopes models
sum((SoutputN[,1]-SoutputN[,2]*1.96)>0) #null
## [1] 1
sum((Soutput1[,1]-Soutput1[,2]*1.96)>0) #smaller effect
## [1] 6
sum((Soutput2[,1]-Soutput2[,2]*1.96)>0) #larger effect
## [1] 3

We now produce different elements of the plot in the main text

#Panel 1, example of raw dataset, points coloured by group
plot(signalsize2~bodysize,
     pch=c(rep(16,n-groupsize),rep(15,groupsize)),
     col=c(rep("grey",n-groupsize),rep("red",groupsize)),
     ylab="Centralised Signal Size",
     xlab="Centralised Body Size",
     cex.axis=1.5,
     cex.lab=1.5,
     las=1)

#Panel 2, Model estimates for random intercept models
par(mar=c(4,5,2,2))

cols=c("#4477AA", "#DDCC77", "#CC6677")
c1<-col2rgb(cols[1])/255
c2<-col2rgb(cols[2])/255

plot(NULL,xlim=c(-0.03,0.05),ylim=c(0,40),xlab="Model estimate for the effect of body size",ylab="Density",cex.lab=1.5,cex.axis=1.5,yaxs="i",las=1)
polygon(density(Ioutput2[,1]),col=rgb(col2rgb(cols[1])[1]/255,col2rgb(cols[1])[2]/255,col2rgb(cols[1])[3]/255,0.1),border=cols[1],lwd=3)
polygon(density(Soutput2[,1]),col=rgb(col2rgb(cols[2])[1]/255,col2rgb(cols[2])[2]/255,col2rgb(cols[2])[3]/255,0.1),border=cols[2],lwd=3)
lines(x=c(-10000,10000),y=c(0,0),lwd=2)

text(-0.018,37,"Random Slopes",cex=1.5,adj=c(0,0.5))
text(-0.018,34,"Random Intercepts",cex=1.5,adj=c(0,0.5))
points(x=-0.023,y=37,pch=22,col=cols[2],bg=rgb(c2[1],c2[2],c2[3],0.2),cex=4)
points(x=-0.023,y=34,pch=22,col=cols[1],bg=rgb(c1[1],c1[2],c1[3],0.2),cex=4)

#Panel 4, Figure of histograms of error distributions in extreme example
par(mar=c(4,5,2,2))
plot(NULL,xlim=c(0,0.04),ylim=c(0,800),xlab="Standard error around the effect of body size",ylab="",cex.lab=1.5,cex.axis=1.5,yaxs="i",las=1)
polygon(density(Ioutput2[,2]),col=rgb(col2rgb(cols[1])[1]/255,col2rgb(cols[1])[2]/255,col2rgb(cols[1])[3]/255,0.1),border=cols[1],lwd=2)
polygon(density(Soutput2[,2]),col=rgb(col2rgb(cols[2])[1]/255,col2rgb(cols[2])[2]/255,col2rgb(cols[2])[3]/255,0.1),border=cols[2],lwd=2)
lines(x=c(-10000,10000),y=c(0,0),lwd=2)

text(0.025,700,"Random Slopes",cex=1.5,adj=c(0,0.5))
text(0.025,650,"Random Intercepts",cex=1.5,adj=c(0,0.5))
points(x=0.022,y=700,pch=22,col=cols[2],bg=rgb(c2[1],c2[2],c2[3],0.2),cex=4)
points(x=0.022,y=650,pch=22,col=cols[1],bg=rgb(c1[1],c1[2],c1[3],0.2),cex=4)

par(xpd=NA)
text(x=-0.007,y=400,labels="Density",srt=90,cex=1.5)

par(xpd=FALSE)