########################################################
#   Resistência de Papel (Montgomery)                  #
##  Split-Plot                                        ##
#    Fatores: método e temperatura                     #
#    Bloco: dia                                        #
########################################################
#
 rm(list=ls(all=TRUE)) #Limpeza do ambiente.
#
#Instale e carregue pacotes necessários para a análise:
#
# install.packages("agricolae"); install.packages("car")
#install.packages("GAD"); install.packages("lattice")
#install.packages("lsmeans");
 require(agricolae); require(car); require(GAD)
 require(lattice); require(lsmeans)
#
#
 e3 <- read.table("split_plot.txt", header = T)
 sapply(e3, class)
 e3
 dim(e3)
#
 e3$bloco <- as.factor(e3$bloco)
 e3$metodo <- as.factor(e3$metodo)
 e3$temperatura <- as.factor(e3$temperatura)
 sapply(e3, class)
#
#
##################################################################
#    Análise Descritiva
#
 with(e3, xyplot(resist ~ temperatura | metodo))
 with(e3, xyplot(resist ~ temperatura | metodo, groups = bloco))
 interaction.plot(e3$temperatura, e3$metodo, e3$resist, type="b", 
 legend=T,trace.label="método", ylab="Resistência do Papel", xlab="temperatura (F)",lwd=2)
#
 boxplot(resist ~ bloco, data=e3,
 ylab="resistência", main="Boxplots dos Dias")
# 
 boxplot(resist ~  metodo, data=e3,
 ylab="resistência", main="Boxplots dos Métodos")
# 
 boxplot(resist ~  bloco*metodo, data=e3,
 ylab="resistência", main="Boxplots dos Blocos e Métodos")
# 
 boxplot(resist ~  temperatura, data=e3,
 ylab="resistência", main="Boxplots das Temperaturas")
# 
 boxplot(resist ~  metodo*temperatura, data=e3,
 ylab="resistência", main="Boxplots dos Temperaturas e Métodos")
#
#####################################################################
#
# Os dois fatores são fixos e o bloco aleatório.
 e3$metodo <- as.fixed(e3$metodo)
 e3$temperatura <- as.fixed(e3$temperatura)
 e3$bloco <- as.random(e3$bloco) 
#
  modelo<-with(e3,sp.plot(block = bloco, pplot = metodo,
       splot = temperatura, Y = resist))
#
 modelo1 <- aov(resist ~ bloco + metodo*temperatura + Error(bloco/metodo),
 data=e3)
 summary(modelo1)
#
# Verificando as suposições do modelo ####
#
# De forma a verificar as suposições do modelo, é necessário usar um termo de erro. 
# Podemos utilizar um modelo não hierárquico para verificar as suposições e também
# para fazer comparações múltiplas.
#    Observe que o teste F tem que ser usado somente no modelo correto.
#
 res.modelo2 <- aov(resist ~ bloco + metodo*temperatura + bloco:metodo, data = e3)
 summary(res.modelo2) # testes F errados
 par(mfrow=c(2,2))
 plot(res.modelo2)
 boxCox(res.modelo2)
#
# Comparações múltiplas
#
#
 tx <- with(e3, interaction(temperatura, metodo))
 amod <- aov(resist ~ tx, data=e3)
 plot(HSD.test(amod, "tx", group=TRUE))
 plot(TukeyHSD(amod, "tx", group=T)
# 
 print(LSD.test(amod, tx, group=TRUE))
 ex01.tu <- TukeyHSD(aov(resist ~ metodo*temperatura, data=e3))
 plot(ex01.tu)
 TukeyHSD(amod, which = "metodo*temperatura")
#
 modelo3<-aov(resist ~ metodo*temperatura, data=e3)
 ex01.tu <- TukeyHSD(modelo3,c("temperatura*metodo"))
 plot(ex01.tu)
#