library(multcomp)

#One-way CRD
y=c(56,48,66,62,83,78,94,93,80,72,83,85)
x=factor(c(rep(1,4),rep(2,4),rep(3,4)))
par(mfrow=c(1,3))
lm1=lm(y~x)
plot(x,y,main="Seedlings",xlab="Insecticide",ylab="# Germinated per 100")
anova(lm1)
qqnorm(lm1$resid,main="Residuals")
qqline(lm1$resid)
r=rnorm(length(y),0,7.18)
qqnorm(r,main="N(0,7.18)")
qqline(r)

# CRBD design
Productivity=c(33,38,39,42,62,35,37,43,47,71,40,42,45,52,74,54,50,55,62,84)
Workshop=factor(c(rep('A',5),rep('B',5),rep('C',5),rep('D',5)))
Attitude=factor(c(rep(1:5,4)))
lmR=lm(Productivity ~ Workshop + Attitude)
anova(lmR)
interaction.plot(Workshop,Attitude,Productivity,ylab="Productivity")


#Mutiple Comparison in one-way layout
SideEffect=c(27,26,21,26,19,13,15,16,15,10,10,11,22,15,21,18,20,18,17,16)
Drug=factor(rep(1:5,each=4))
lmod=aov(SideEffect~Drug)
anova(lmod)
contr=rbind("l1"=c(1,-1/4,-1/4,-1/4,-1/4),"l2"=c(0,1,-1,0,0),"l3"=c(0,0,0,1,-1),l4=c(0,1/2,1/2,-1/2,-1/2))
mc=glht(lmod,linfct=mcp(Drug=contr))
summary(mc,test=adjusted(type="bonferroni"))
confint(mc,calpha=qt(0.975,15)) #Fisher's Method
confint(mc,calpha=qt(1-0.05/8,15)) #Fishers Bonferroni's Method
confint(mc,calpha=sqrt(4*qf(1-0.05,4,15))) #Scheffe's Method
par(mfrow=c(1,3))
plot(confint(mc,calpha=qt(0.975,15)),main="Fisher",xlab="95% Confidence Interval")
plot(confint(mc,calpha=qt(1-0.05/8,15)),main="Bonferroni",xlab="95% Confidence Interval")
plot(confint(mc,calpha=sqrt(4*qf(1-0.05,4,15))),main="Scheffe",xlab="95% Confidence Interval")

#Dunnett's Procedure
n=tapply(Drug,Drug,length)
Dunn=glht(lmod,linfct=mcp(Drug="Dunnett"),alternative="less")
summary(Dunn)
confint(Dunn)
par(mfrow=c(1,1))
plot(confint(Dunn))

# Multiple Pairwise Comparisons
mupc=glht(lmod,linfct=mcp(Drug="Tukey"))
summary(mupc)
summary(mupc,test=adjusted("bonferroni"))

confint(mupc,calpha=qt(1-0.05/2,15)) #Fisher's Method
confint(mupc,calpha=qt(1-0.05/20,15)) #Bonferroni's Method
confint(mupc,calpha=sqrt(4*qf(1-0.05,4,15))) #Scheffe's Method)
confint(mupc) #Tukey's Method
par(mfrow=c(1,4))
plot(confint(mupc,calpha=qt(1-0.05/2,15)))
plot(confint(mupc,calpha=qt(1-0.05/20,15)))
plot(confint(mupc,calpha=sqrt(4*qf(1-0.05,4,15))))
plot(confint(mupc))
Tukey.cld=cld(mupc)
par(mfrow=c(1,1))
old.par=par(mai=c(1,1,1.25,1))
plot(Tukey.cld)
par(old.par)

# Tukey's Multiple Comparison for RCBD
Productivity=c(33,38,39,42,62,35,37,43,47,71,40,42,45,52,74,54,50,55,62,84)
Workshop=factor(c(rep('A',5),rep('B',5),rep('C',5),rep('D',5)))
Attitude=factor(c(rep(1:5,4)))
lmR=lm(Productivity ~ Workshop + Attitude)
mcR=glht(lmR,linfct=mcp(Workshop="Tukey"))
summary(mcR)
mcR.cld=cld(mcR)
par(mfrow=c(1,1))
old.par=par(mai=c(1,1,1.25,1))
plot(mcR.cld)
par(old.par)
mcRW=glht(lmR,linfct=mcp(Workshop="Williams"))
plot(confint(mcRW,calpha=sqrt(3*qf(1-0.05,3,16))))

#Tukey's Multiple Comparison in Two-way Anova with insignificant interaction
Yield=c(49,39,50,55,43,38,53,48,55,41,65,58,53,42,85,73,66,68,85,92,69,62,85,99)
Variety=factor(rep(1:3,each=8))
Pesticide=factor(rep(rep(1:4,each=2),times=3))
interaction.plot(Variety,Pesticide,Yield)
amod=aov(Yield~ Variety * Pesticide)
anova(amod)
mcV=glht(amod,linfct=mcp(Variety="Tukey"))
summary(mcV)
(mcV.cld=cld(mcV))
mcP=glht(amod,linfct=mcp(Pesticide="Tukey"))
summary(mcP)
(mcP.cld=cld(mcP))

#Tukey's Multiple Comparisons in Two-Factor ANOVA with significant Interaction 
Data=read.csv(file="Orange.csv")
Data$PH=factor(Data$PH)
Data$calcium=factor(Data$calcium)
par(mfrow=c(1,2))
interaction.plot(Data$PH,Data$calcium,Data$Tr_Diam,fun=mean)
interaction.plot(Data$calcium,Data$PH,Data$Tr_Diam,fun=mean)
lm1=lm(Tr_Diam~PH+calcium + PH:calcium,data=Data)
mcPH=glht(lm1,linfct=mcp(PH="Tukey"))
summary(mcPH)
mccalcium=glht(lm1,linfct=mcp(calcium="Tukey"))
summary(mccalcium)