pdf("pbahasse.pdf") pbahasse <- read.fwf(file="http://www.stat.ufl.edu/~winner/computing/sas/pbahasse.dat", width=c(8,8,8,8,8,8), col.names=c("score","tourney","rep","bowler","oil","touney_oil")) attach(pbahasse) tourney <- factor(tourney, levels=1:15) bowler <- factor(bowler,levels=1:37) oil <- factor(oil, levels=1:5) score <- score/100 pba.aov <- aov(score ~ oil*bowler*tourney) summary(pba.aov) library(lmerTest) pba.mod1 <- lmer(score ~ oil + (1|bowler) + (1|oil:bowler) + (1|oil:tourney) + (1|bowler:oil:tourney)) summary(pba.mod1) anova(pba.mod1) if(require(pbkrtest)) anova(pba.mod1, ddf = "Kenward-Roger") lsmeans(pba.mod1) difflsmeans(pba.mod1)