Comparing R TukeyHSD results with Statistics, William L. Hays p 581-583

48872315.R

# StackOverflow question: 
# https://stackoverflow.com/questions/48872315/getting-wrong-p-values-for-tukey-test-for-one-way-mixed-effect-anova/

subject=c(1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 5, 6, 6, 6, 
          6, 7, 7, 7, 7, 8, 8, 8, 8, 9, 9, 9, 9, 10, 10, 10, 10)

treatment=c("a1", "a2", "a3", "a4", "a1", "a2", "a3", "a4", "a1", "a2", "a3", 
            "a4", "a1", "a2", "a3", "a4", "a1", "a2", "a3", "a4", "a1", "a2", 
            "a3", "a4", "a1", "a2", "a3", "a4", "a1", "a2", "a3", "a4", "a1", 
            "a2", "a3", "a4", "a1", "a2", "a3", "a4")

response=c(11, 9, 5, 17, 14, 12, 10, 18, 15, 7, 6, 21, 17, 10, 13, 22, 15, 7, 6, 
           15, 14, 8, 13, 22, 9, 6, 7, 15, 17, 11, 10, 19, 10, 13, 14, 23, 12, 
           8, 11, 20)

dataFrame=data.frame(subject, treatment, response)

library(nlme)
library(multcomp)
set.seed(1)
model <- lme(response~treatment, random=~1|subject, data=dataFrame, method="ML")
anova.lme(model)

glht.out <- glht(model, mcp(treatment="Tukey"), alternative="two.sided")
summary(glht.out)

?summary.glht
summary(glht.out, test = adjusted(type = "single-step"))
summary(glht.out, test = adjusted(type = "none"))

aovmodel <- aov(model)
summary(aovmodel)
TukeyHSD(aovmodel, "treatment")
# https://github.com/SurajGupta/r-source/blob/master/src/library/stats/R/TukeyHSD.R#L72
MSE <- 8.15
center <- -3.9/sqrt(MSE/10)
ptukey(abs(center), 4, 36, lower.tail=FALSE)