library(faraway) data(pulp) op <- options(contrasts=c("contr.sum", "contr.poly")) lmod <- aov(bright ~ operator, pulp) summary(lmod) coef(lmod) options(op) (0.447-0.106)/5 library(lme4) mmod <- lmer(bright ~ 1+(1|operator), pulp) summary(mmod) smod <- lmer(bright ~ 1+(1|operator), pulp, method="ML") summary(smod) nullmod <- lm(bright ~ 1, pulp) as.numeric(2*(logLik(smod)-logLik(nullmod))) pchisq(2.5684,1,lower=FALSE) y <- simulate(nullmod) lrstat <- numeric(1000) for(i in 1:1000){ y <- unlist(simulate(nullmod)) bnull <- lm(y ~ 1) balt <- lmer(y ~ 1 + (1|operator),pulp,method="ML") lrstat[i] <- as.numeric(2*(logLik(balt)-logLik(bnull))) } mean(lrstat < 0.00001) mean(lrstat > 2.5684) sqrt(0.02*0.98/1000) ranef(mmod)$operator (cc <- model.tables(lmod)) cc[[1]]$operator/ranef(mmod)$operator fixef(mmod)+ranef(mmod)$operator qqnorm(resid(mmod),main="") plot(fitted(mmod),resid(mmod),xlab="Fitted",ylab="Residuals") abline(0,0) data(penicillin) summary(penicillin) op <- options(contrasts=c("contr.sum", "contr.poly")) lmod <- aov(yield ~ blend + treat, penicillin) summary(lmod) coef(lmod) mmod <- lmer(yield ~ treat + (1|blend), penicillin) summary(mmod) options(op) ranef(mmod)$blend anova(mmod) amod <- lmer(yield ~ treat + (1|blend), penicillin, method="ML") nmod <- lmer(yield ~ 1 + (1|blend), penicillin, method="ML") anova(amod,nmod) lrstat <- numeric(1000) for(i in 1:1000){ ryield <- unlist(simulate(nmod)) nmodr <- lmer(ryield ~ 1 + (1|blend), penicillin, method="ML") amodr <- lmer(ryield ~ treat + (1|blend), penicillin, method="ML") lrstat[i] <- 2*(logLik(amodr)-logLik(nmodr)) } plot(qchisq((1:1000)/1001,3),sort(lrstat),xlab=expression(chi[3]^2),ylab="Simulated LRT") abline(0,1) mean(lrstat > 4.05) rmod <- lmer(yield ~ treat + (1|blend), penicillin) nlmod <- lm(yield ~ treat, penicillin) 2*(logLik(rmod)-logLik(nlmod,REML=TRUE)) lrstatf <- numeric(1000) for(i in 1:1000){ ryield <- unlist(simulate(nlmod)) nlmodr <- lm(ryield ~ treat, penicillin) rmodr <- lmer(ryield ~ treat + (1|blend), penicillin) lrstatf[i] <- 2*(logLik(rmodr)-logLik(nlmodr,REML=TRUE)) } mean(lrstatf < 0.00001) cs <- lrstatf[lrstatf > 0.00001] ncs <- length(cs) plot(qchisq((1:ncs)/(ncs+1),1),sort(cs),xlab=expression(chi[1]^2),ylab="Simulated LRT") abline(0,1) mean(lrstatf > 2.7629) data(irrigation) summary(irrigation) lmod <- lmer(yield ~ irrigation * variety + (1|field) +(1|field:variety),data=irrigation) logLik(lmod) lmodr <- lmer(yield ~ irrigation * variety + (1|field),data=irrigation) logLik(lmodr) summary(lmodr) anova(lmodr) plot(fitted(lmodr),resid(lmodr),xlab="Fitted",ylab="Residuals") qqnorm(resid(lmodr),main="") mod <- lm(yield ~ irrigation * variety, data=irrigation) anova(mod) data(eggs) summary(eggs) cmod <- lmer(Fat ~ 1 + (1|Lab) + (1|Lab:Technician) + (1|Lab:Technician:Sample), data=eggs) summary(cmod) cmodr <- lmer(Fat ~ 1 + (1|Lab) + (1|Lab:Technician), data=eggs) anova(cmod,cmodr) VarCorr(cmodr) lrstat <- numeric(1000) for(i in 1:1000){ rFat <- unlist(simulate(cmodr)) nmod <- lmer(rFat ~ 1 + (1|Lab) + (1|Lab:Technician), data=eggs) amod <- lmer(rFat ~ 1 + (1|Lab) + (1|Lab:Technician) + (1|Lab:Technician:Sample), data=eggs) lrstat[i] <- 2*(logLik(amod)-logLik(nmod)) } mean(lrstat < 0.00001) 2*(logLik(cmod)-logLik(cmodr)) mean(lrstat > 1.6034) data(abrasion) matrix(abrasion$material,4,4) lmod <- aov(wear ~ material + run + position, abrasion) summary(lmod) mmod <- lmer(wear ~ material + (1|run) + (1|position), abrasion) anova(mmod) summary(mmod) data(jsp) jspr <- jsp[jsp$year==2,] plot(jitter(math)~jitter(raven),data=jspr,xlab="Raven score",ylab="Math score") boxplot(math ~ social, data=jspr,xlab="Social class",ylab="Math score") glin <- lm(math ~ raven*gender*social,jspr) anova(glin) glin <- lm(math ~ raven*social,jspr) anova(glin) glin <- lm(math ~ raven+social,jspr) summary(glin) table(jspr$school) mmod <- lmer(math ~ raven*social*gender+(1|school)+(1|school:class),data=jspr) anova(mmod) jspr$craven <- jspr$raven-mean(jspr$raven) mmod <- lmer(math ~ craven*social+(1|school)+(1|school:class),jspr) summary(mmod) qqnorm(resid(mmod),main="") plot(fitted(mmod),resid(mmod),xlab="Fitted",ylab="Residuals") qqnorm(ranef(mmod)$school[,1],main="School effects") qqnorm(ranef(mmod)$"school:class"[,1],main="Class effects") adjscores <- ranef(mmod)$school[,1] rawscores <- coef(lm(math ~ school-1,jspr)) rawscores <- rawscores-mean(rawscores) plot(rawscores,adjscores) sint <- c(9,14,29) text(rawscores[sint],adjscores[sint]+0.2,c("9","15","30")) schraven <- lm(raven ~ school, jspr)$fit mmodc <- lmer(math ~ craven*social+schraven*social+(1|school)+(1|school:class),jspr) anova(mmodc)