Peril 3: Number of levels of the random effect

Load required R packages

library(lme4)

Define parameters used in the simulations

#Number of grouping levels
groups<-c(seq(2,20,1),22,24,26,28,30,35,40)
#Total sample size
n<-200
#Standard deviation of the random effects distribution
gr.dev<-3
#Population intercept
intercept<-50
#Standard deviation for the residual deviance
resid.dev<-3
#Effect of the covariate
x.eff<-1

#Values of the explanatory variable
x<-round(runif(n,min=-20,max=20))

#Number of simulations to be conducted
reps<-100

Now we create empty matrices to store the results

rvar<-matrix(NA,nr=reps,nc=length(groups))
int_est<-matrix(NA,nr=reps,nc=length(groups))
x_est<-matrix(NA,nr=reps,nc=length(groups))
int_err<-matrix(NA,nr=reps,nc=length(groups))
x_err<-matrix(NA,nr=reps,nc=length(groups))

Now we can run 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.

#loop over number of simulations
for(r in 1:reps){
  
 #loop over the number of random effect levels  
 for(i in 1:length(groups)){
  
  #simulate group data
  group<-sample(1:groups[i],200,replace=TRUE)
  group.means<-rnorm(groups[i],mean=0,sd=gr.dev)
  
  #Simulate response variable
  y<-numeric()
  for(j in 1:n){
    y[j]<-rnorm(1,intercept+x.eff*x[j]+group.means[group[j]],resid.dev)
  }
  
  #Fit mixed model in lme4
  mod<-lmer(y~x+(1|group))
  
  #Quantify model estimates that we want to keep
  rvar[r,i]<-sqrt(unlist(VarCorr(mod)[1]))
  int_est[r,i]<-fixef(mod)[1]
  x_est[r,i]<-fixef(mod)[2]
  int_err[r,i]<-summary(mod)$coefficients[1,2]
  x_err[r,i]<-summary(mod)$coefficients[2,2]
  
 } #end groups loop
} #end reps loop

Now we can plot the results, as used in the main text

par(mfrow=c(2,1),mar=c(4,6,2,2))
boxplot(rvar,col="light grey",range=0,las=1,lty=1,names=groups)
lines(x=c(-20,100),y=c(3,3),lwd=4,col="red")
mtext("Random Effect Variance (SD)", side=2,line=3.4,cex=1.5)
mtext("a)",side=3,line=-1.2,adj=0.003,cex=1.5)
mtext("Number of levels of random effect",side=1,line=2.5,cex=1.5)
boxplot(x_err,col="light grey",range=0,las=1,lty=1,names=groups)
mtext("Standard Error", side=2,line=3.4,cex=1.5)
mtext("b)",side=3,line=-1.2,adj=0.003,cex=1.5)
mtext("Number of levels of random effect",side=1,line=2.5,cex=1.5)