Peril 1: (G)LMM p-values can be anti-conservative

Load required R packages

library(nlme)
library(lme4)
library(lmerTest)

PART 1: Low replication ***

Set up empty vectors to store results

cis<-numeric(1000)
ltests<-numeric(1000)
AICcheck.ML<-numeric(1000)
cis.lmer<-numeric(1000)
ltests.lmer<-numeric(1000)
AICcheck.lmer.ML<-numeric(1000)

Define parameters

#Define number of clones
clones<-6

#Number of samples per clone
perclone<-2

Run 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 for 1000 iterations
for(i in 1:1000){

 #generate clone factor
 clone<-gl(clones,perclone)
 #generate treatment factor
 treat<-rep(gl(2,perclone/2),clones)
 #create group-level means for response variable
 y.mean<-rep(rnorm(clones,0,1),each=perclone)
 #create values for response variable
 y.obs<-rnorm(perclone*clones,y.mean,1)
 #treatment model
 m1<-lme(y.obs~treat,random=~1|clone, control = lmeControl(opt = "optim"))
 #null model
 m2<-lme(y.obs~1,random=~1|clone, control = lmeControl(opt = "optim"))
 #confidence interval overlaps zero or not
 t_ci<-abs(summary(m1)$tTable[2,1])-abs(summary(m1)$tTable[2,2]*1.96)
 cis[i]<-as.numeric(t_ci>0)
 #fits same model and then does anova with different "control" function
 m3<-lme(y.obs~treat,random=~1|clone,method="ML", control = lmeControl(opt = "optim"))
 m4<-update(m3,~.-treat, control = lmeControl(opt = "optim"))
 #p value from model anova
 ltests[i]<-anova(m3,m4)$"p-value"[2]
 #AICs
 AICcheck.ML[i]<-AIC(m3)-AIC(m4)
 #lme4 version
 m1.lmer<-lme4::lmer(y.obs~treat+(1|clone))
 m2.lmer<-lme4::lmer(y.obs~1+(1|clone))
 #confidence interval overlaps zero or not from lme4
 t_ci<-abs(summary(m1.lmer)$coefficients[2,1])-abs(summary(m1.lmer)$coefficients[2,2]*1.96)
 cis.lmer[i]<-as.numeric(t_ci>0)
 #Refit model with REML=F
 m3.lmer<-lme4::lmer(y.obs~treat+(1|clone),REML=F)
 m4.lmer<-lme4::lmer(y.obs~1+(1|clone),REML=F)
 #AICs
 AICcheck.lmer.ML[i]<-AIC(m3.lmer)-AIC(m4.lmer)
 #anova of treatment and null model
 ltests.lmer[i]<-anova(m3.lmer,m4.lmer)$"Pr(>Chisq)"[2]
}

Now we work out the Type I error rate and print them in the order: - whether 95% confidence intervals around treatment effect overlap zero (nlme) - likelihood ratio test (nlme) - AIC (nlme) - whether 95% confidence intervals around treatment effect overlap zero (lme4) - likelihood ratio test (lme4) - AIC (lme4)

Alpha for likelihood ratio test is 0.05 Threshold used for AIC is 2

A<-mean(cis)
B<-length(ltests[ltests>0&ltests<0.05])/length(ltests[ltests>0])
C<-length(AICcheck.ML[ltests>0&AICcheck.ML<(-2)])/length(AICcheck.ML[ltests>0])
D<-mean(cis.lmer)
E<-length(ltests.lmer[ltests.lmer>0&ltests.lmer<0.05])/length(ltests.lmer[ltests.lmer>0])
F<-length(AICcheck.lmer.ML[ltests.lmer>0&AICcheck.lmer.ML<(-2)])/length(AICcheck.lmer.ML[ltests.lmer>0])

RES1<-list(A,B,C,D,E,F)

print(RES1)
## [[1]]
## [1] 0.122
## 
## [[2]]
## [1] 0.105
## 
## [[3]]
## [1] 0.1
## 
## [[4]]
## [1] 0.122
## 
## [[5]]
## [1] 0.105
## 
## [[6]]
## [1] 0.1

PART 2: High replication replication ***

Set up empty vectors to store results

cis<-numeric(1000)
ltests<-numeric(1000)
AICcheck.ML<-numeric(1000)
cis.lmer<-numeric(1000)
ltests.lmer<-numeric(1000)
AICcheck.lmer.ML<-numeric(1000)

Define parameters - note the sample size for each clone is now 10

#Define number of clones
clones<-6

#Number of samples per clone
perclone<-10

Run 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 for 1000 iterations
for(i in 1:1000){
 #generate clone factor
 clone<-gl(clones,perclone)
 #generate treatment factor
 treat<-rep(gl(2,perclone/2),clones)
 #create group-level means for response variable
 y.mean<-rep(rnorm(clones,0,1),each=perclone)
 #create values for response variable
 y.obs<-rnorm(perclone*clones,y.mean,1)
 #treatment model
 m1<-lme(y.obs~treat,random=~1|clone, control = lmeControl(opt = "optim"))
 #null model
 m2<-lme(y.obs~1,random=~1|clone, control = lmeControl(opt = "optim"))
 #confidence interval overlaps zero or not
 t_ci<-abs(summary(m1)$tTable[2,1])-abs(summary(m1)$tTable[2,2]*1.96)
 cis[i]<-as.numeric(t_ci>0)
 #fits same model and then does anova with different "control" function
 m3<-lme(y.obs~treat,random=~1|clone,method="ML", control = lmeControl(opt = "optim"))
 m4<-update(m3,~.-treat, control = lmeControl(opt = "optim"))
 #p value from model anova
 ltests[i]<-anova(m3,m4)$"p-value"[2]
 #AICs
 AICcheck.ML[i]<-AIC(m3)-AIC(m4)
 #lme4 version
 m1.lmer<-lme4::lmer(y.obs~treat+(1|clone))
 m2.lmer<-lme4::lmer(y.obs~1+(1|clone))
 #confidence interval overlaps zero or not from lme4
 t_ci<-abs(summary(m1.lmer)$coefficients[2,1])-abs(summary(m1.lmer)$coefficients[2,2]*1.96)
 cis.lmer[i]<-as.numeric(t_ci>0)
 #Refit model with REML=F
 m3.lmer<-lme4::lmer(y.obs~treat+(1|clone),REML=F)
 m4.lmer<-lme4::lmer(y.obs~1+(1|clone),REML=F)
 #AICs
 AICcheck.lmer.ML[i]<-AIC(m3.lmer)-AIC(m4.lmer)
 #anova of treatment and null model
 ltests.lmer[i]<-anova(m3.lmer,m4.lmer)$"Pr(>Chisq)"[2]
}

Now we work out the Type I error rate and print them in the order: - whether 95% confidence intervals around treatment effect overlap zero (nlme) - likelihood ratio test (nlme) - AIC (nlme) - whether 95% confidence intervals around treatment effect overlap zero (lme4) - likelihood ratio test (lme4) - AIC (lme4)

Alpha for likelihood ratio test is 0.05 Threshold used for AIC is 2

A2<-mean(cis)
B2<-length(ltests[ltests>0&ltests<0.05])/length(ltests[ltests>0])
C2<-length(AICcheck.ML[ltests>0&AICcheck.ML<(-2)])/length(AICcheck.ML[ltests>0])
D2<-mean(cis.lmer)
E2<-length(ltests.lmer[ltests.lmer>0&ltests.lmer<0.05])/length(ltests.lmer[ltests.lmer>0])
F2<-length(AICcheck.lmer.ML[ltests.lmer>0&AICcheck.lmer.ML<(-2)])/length(AICcheck.lmer.ML[ltests.lmer>0])

RES2<-list(A2,B2,C2,D2,E2,F2)

print(RES2)
## [[1]]
## [1] 0.059
## 
## [[2]]
## [1] 0.057
## 
## [[3]]
## [1] 0.051
## 
## [[4]]
## [1] 0.059
## 
## [[5]]
## [1] 0.057
## 
## [[6]]
## [1] 0.033