#install.packages("WRS2")
library(WRS2)
data("goggles")
attach(goggles)
levels(alcohol)


contrasts(alcohol)="contr.treatment"
mod0=lm(attractiveness~alcohol,goggles)
coef(mod0)
contrasts(alcohol)

aggregate(attractiveness~alcohol,goggles,FUN=mean)

contrasts(alcohol)="contr.sum"
mod0=lm(attractiveness~alcohol)
coef(mod0)

boxplot(attractiveness~alcohol,goggles)
mod0=lm(attractiveness~alcohol,goggles)
anova0=anova(mod0)
anova(mod0)

mod1=aov(attractiveness~alcohol,goggles)
anova1=summary(mod1)

summary(mod1)

means=aggregate(attractiveness~alcohol,goggles,FUN=mean)

n=nrow(goggles)
mse=anova0$`Mean Sq`[2]
Dnone2=means$attractiveness[2]-means$attractiveness[1]
Dnone4=means$attractiveness[3]-means$attractiveness[1]
D24=means$attractiveness[3]-means$attractiveness[2]

q=qtukey(0.05,3,48-3,lower.tail = F)

ni=aggregate(attractiveness~alcohol,goggles,length)
nnone=ni$attractiveness[1]
n2=ni$attractiveness[2]
n4=ni$attractiveness[3]

#1
Dnone2+c(-1,1)*1/sqrt(2)*q*sqrt(mse*(1/nnone+1/n2))

#2
Dnone4+c(-1,1)*1/sqrt(2)*q*sqrt(mse*(1/nnone+1/n4))

#3
D24+c(-1,1)*1/sqrt(2)*q*sqrt(mse*(1/n2+1/n4))