#
#   Exemplo: Bateria - Dois Fatores
#            Fatores: tipo (3) e temperatura (3)
#            4 replicas
#------------------------------------------------------------------------
 rm(list=ls(all=TRUE))
 require(multcomp)
 require(agricolae)
 require(nortest)
#require(car)
#
 dados=read.table("dados2_doisFatores_bateria.txt", header=TRUE)
 names(dados)
 sapply(dados,class)
#
 dados$tipo=as.factor(dados$tipo)
 dados$temperatura=as.factor(dados$temperatura)
#
 interaction.plot(dados$temperatura, dados$tipo, dados$y, type="b", ylab="tempo de vida da bateria", lwd=2)
 interaction.plot(dados$temperatura, dados$tipo, dados$y, type="b", ylab="tempo de vida da bateria", lwd=2)
 boxplot(y ~ temperatura * tipo, data=dados,
 ylab="Tempo de vida", main="Boxplots Temperatura e Tipo")
 boxplot(y ~ temperatura , data=dados,
 ylab="Tempo de vida", main="Boxplots Temperatura")
 boxplot(y ~ tipo, data=dados,
 ylab="Tempo de vida", main="Boxplots Tipo")
#
 modelo=aov(y~tipo*temperatura, data=dados)
 summary(modelo)
#    
################################
#    Análise dos resíduos
#
 y.hat=modelo$fitted.values
 e=modelo$residuals
 par(mfrow=c(2,2))
 plot(modelo)
#
 tipo=as.numeric(dados[,2])
 temperatura=as.numeric(dados[,3])
# qqnorm(e)
# qqline(e, col="red")
 hist(e)
#
# Testes de normalidade:
 shapiro.test(e)
 ad.test(e)
#
# Testes de homogeneidade da variância
 bartlett.test(y~temperatura * tipo, data=dados)
#
################################
# Comparações m?ltiplas
#
 model.tables(modelo, ty="means")
 library(agricolae)
 tx <- with(dados, interaction(temperatura, tipo))
 amod <- aov(y ~ tx, data=dados)
 plot(HSD.test(amod, "tx", group=TRUE))
#
# print(LSD.test(modelo, "tipo", group=FALSE))
# print(LSD.test(modelo, "temperatura", group=FALSE))
 print(LSD.test(modelo, c("temperatura","tipo"), group=TRUE))
 ex01.tu <- TukeyHSD(modelo,c("temperatura","tipo","temperatura*tipo"))
 plot(ex01.tu)
#
###################################################################################
#
#   Exemplo: Impureza - Dois Fatores sem replicação
#            Fatores: temperatura (3) e pressão (5)
#            
#------------------------------------------------------------------------
 rm(list=ls(all=TRUE))
 require(multcomp)
 require(agricolae)
 require(nortest)
#require(car)
#
 dados1=read.table("impureza.txt", header=TRUE)
 names(dados1)
 sapply(dados1,class)
#
 dados1$Temp=as.factor(dados1$Temp)
 dados1$pressao=as.factor(dados1$pressao)
#
 interaction.plot(dados1$Temp, dados1$pressao, dados1$Impureza, type="b", 
                  legend=T,trace.label="pressão", ylab="% impureza",xlab="temperatura (F)", lwd=2)
#
 tukeys.add.test(log(dados1$Impureza),dados1$Temp,dados1$pressao)
#
 boxplot(Impureza ~ Temp , data=dados1,
 ylab="% Impureza", main="Boxplots Temperatura")
 boxplot(Impureza ~ pressao , data=dados1,
 ylab="% Impureza", main="Boxplots Pressão")
#
 modelo1<-aov(dados1$Impureza ~ dados1$Temp + dados1$pressao)
 summary(modelo1)
################################
#    Análise dos resíduos
#
 y.hat=modelo1$fitted.values
 e=modelo1$residuals
 par(mfrow=c(2,2))
 plot(modelo1)
#
# Testes de normalidade:
 shapiro.test(e)
 ad.test(e)
# Testes de homogeneidade da variância ---> não funciona pois não temos replicas.
# bartlett.test(split(as.numeric(dados1$Impureza), list(dados1$Temp, dados1$pressao)))
#
###################################################################
# 
 colnames(dados1) <- c("Impureza","Temperatura", "Pressao")
 modelo2<-aov(log(Impureza) ~ Temperatura + Pressao, data=dados1)
 summary(modelo2)
################################
#    Análise dos resíduos
#
 y.hat=modelo2$fitted.values
 e=modelo2$residuals
 par(mfrow=c(2,2))
 plot(modelo2)
#
# Testes de normalidade:
 shapiro.test(e)
###########################################################################
#  Comparações Mútliplas
#
 ex01.tu <- TukeyHSD(modelo2,dados1$Temperatura)
 plot(ex01.tu)
 ex02.tu <- TukeyHSD(modelo2,dados1$Pressao)
 plot(ex02.tu)
#
#
################################################################################
# Tukey's Test for Additivity:  
##
# The following is an R function to perform Tukey's test of additivity.   
# This function was originally written by Jim Robison-Cox.
# Just copy it into R and run the function with the name of   
# your response and your factors A and B, as shown in the  
# example below. 
# source("C:/Users/Mihinda/Desktop/tukey.txt ")
# data=read.table("C:/Users/Mihinda/Desktop/ex515.txt", header=T) #the data file
# tukeys.add.test(y = data$Data, A = data$Row, B = data$Column)
#                                                   
 tukeys.add.test <- function(y,A,B){
## Y is the response vector
## A and B are factors used to predict the mean of y
## Note the ORDER of arguments: Y first, then A and B
   dname <- paste(deparse(substitute(A)), "and", deparse(substitute(B)),
                  "on",deparse(substitute(y)) )
   A <- factor(A); B <- factor(B)
   ybar.. <- mean(y)
   ybari. <- tapply(y,A,mean)
   ybar.j <- tapply(y,B,mean)
   len.means <- c(length(levels(A)), length(levels(B)))
   SSAB <- sum( rep(ybari. - ybar.., len.means[2]) * 
                rep(ybar.j - ybar.., rep(len.means[1], len.means[2])) *
                tapply(y, interaction(A,B), mean))^2 / 
                  ( sum((ybari. - ybar..)^2) * sum((ybar.j - ybar..)^2))
   aovm <- anova(lm(y ~ A+B))
   SSrem <- aovm[3,2] - SSAB
   dfdenom <- aovm[3,1] - 1
    STATISTIC <- SSAB/SSrem*dfdenom
    names(STATISTIC) <- "F"
    PARAMETER <- c(1, dfdenom)
    names(PARAMETER) <- c("num df", "denom df")
    D <- sqrt(SSAB/  ( sum((ybari. - ybar..)^2) * sum((ybar.j - ybar..)^2)))
    names(D) <- "D estimate"
    RVAL <- list(statistic = STATISTIC, parameter = PARAMETER, 
        p.value = 1 - pf(STATISTIC, 1,dfdenom), estimate = D,
        method = "Tukey's one df F test for Additivity", 
        data.name = dname)
    attr(RVAL, "class") <- "htest"
    return(RVAL)
   } 
##########################################################################
#   Exemplo: Engarrafamento - Três Fatores
#            Fatores: carbonatação (3), Pressão (2) e Velocidade (2)
#            2 replicas
#------------------------------------------------------------------------
 rm(list=ls(all=TRUE))
 require(multcomp)
 require(agricolae)
 require(nortest)
#require(car)
#
 dados2=read.table("garrafa.txt", header=TRUE)
 names(dados2)
 sapply(dados2,class)
 dados2
#
 dados2$carbona=as.factor(dados2$carbona)
 dados2$pressao=as.factor(dados2$pressao)
 dados2$velocida=as.factor(dados2$velocida)
#
 interaction.plot(dados2$carbona, dados2$pressao, dados2$resp, type="b", ylab="desvio nível (mm)", 
 trace.label="pressão", xlab="carbonatação",lwd=2)
 interaction.plot(dados2$carbona, dados2$velocida, dados2$resp, type="b", ylab="desvio nível (mm)", 
 trace.label="velocidade", xlab="pressão",lwd=2)
 interaction.plot(dados2$pressao, dados2$velocida, dados2$resp, type="b", ylab="desvio nível (mm)", 
 trace.label="velocidade", xlab="pressão",lwd=2)
#
 boxplot(resp ~ carbona * velocida * pressao, data=dados2,
 ylab="desvio nível (mm)", main="Boxplots")
 abline(h=0,col=2)
 boxplot(resp ~ carbona * pressao, data=dados2,
 ylab="desvio nível (mm)", main="Boxplots")
 abline(h=0,col=2)
#
 boxplot(resp ~ carbona , data=dados2,
 ylab="desvio nível (mm)", main="Boxplots Carbonatação")
 abline(h=0,col=2)
 boxplot(resp ~ velocida , data=dados2,
 ylab="desvio nível (mm)", main="Boxplots Velocidade")
 abline(h=0,col=2)
 boxplot(resp ~ pressao , data=dados2,
 ylab="desvio nível (mm)", main="Boxplots Pressao")
 abline(h=0,col=2)
#
 modelo1=aov(resp~carbona * velocida * pressao, data=dados2)
 summary(modelo1)
 modelo2=aov(resp~carbona * pressao + velocida * carbona + velocida:pressao, data=dados2)
 summary(modelo2)
#    
################################
#    Análise dos resíduos
#
 y.hat=modelo2$fitted.values
 e=modelo2$residuals
 par(mfrow=c(2,2))
 plot(modelo1)
#
 hist(e)
#
# Testes de normalidade:
 shapiro.test(e)
 ad.test(e)
#
# Testes de homogeneidade da variância
# bartlett.test(y~temperatura * tipo, data=dados)
#
################################
# Comparações m?ltiplas
#
 model.tables(modelo2, ty="means")
 library(agricolae)
 tx <- with(dados2, interaction(carbona, pressao))
 amod <- aov(resp ~ tx, data=dados2)
 plot(HSD.test(amod, "tx", group=TRUE), xlab="carb-pressâo",ylab="desvio nível (mm)")
 lines( abline(h=0, col=12))
#
 tx1 <- with(dados2, velocida)
 amod <- aov(resp ~ tx1, data=dados2)
 plot(HSD.test(amod, "tx1", group=TRUE), xlab="velocidade",ylab="desvio nível (mm)")
#
# print(LSD.test(modelo, "tipo", group=FALSE))
# print(LSD.test(modelo, "temperatura", group=FALSE))
 print(LSD.test(modelo2, c("pressao","carbona"), group=TRUE))
 ex01.tu <- TukeyHSD(modelo2,c("pressao","carbona","pressao*carbona"))
 plot(ex01.tu)
#############################################################