#' ---
#' title: "LCE0602 - Estatística Experimental - Aula03"
#' author: ""
#' date: "08/2023"
#' ---
#' 

#' # Importando um arquivo csv
#' 
#' ## Entrada dos dados 
#' Copiar o arquivo aula2.csv para seu computador
dados <- read.csv2("aula2.csv")
dados$trat<-as.factor(dados$trat)

summary(dados)

#' ## Análise exploratória
library(ggplot2)

ggplot(dados, aes(x = trat, y = y)) +
  geom_point()


#' # Análise de variância  
modelo = aov(y ~ trat, dados)
anova(modelo)


#' # Análise de variância: pressuposições do modelo estatístico
#' Aditividade
#' Independência
#' Hocedasticia
#' Normalidade   

#' ## Resíduos simples
(res <- residuals(modelo))

#' ## Resíduos studentizados
(res_Stud <- rstandard(modelo))

#' Resíduos
round(head(data.frame(res,res_Stud)),5)

#' # Verificação das pressuposições
#' 
#' ## Homogeneidade de variâncias
boxplot(res_Stud)

ggplot(dados, aes(x = trat, y = res_Stud)) +
  geom_point() +
  xlab("Genótipos de milho") +
  ylab("resíduos studentizados")

#' Há uma aparente homogeneidade de variâncias.
#'
#'
#' ### Teste de hipóteses
#' $H_0:$ Há homogeneidade de variâncias *versus* 
#' $H_1:$ Não há homogeneidade de variâncias

library(lawstat)
with(dados,
     levene.test(y, trat, location = "mean"))

#' Como o valor-p é maior do que $\alpha=0,05$, não há evidências
#' para rejeitarmos $H_0$, considerando-se o nível de 5\% de 
#' significância.
#'
#'
#' ## Normalidade dos resíduos

qqnorm(res_Stud)
qqline(res_Stud, col=2)

#' ### Teste de hipóteses
#' $H_0:$ Os erros seguem uma distribuição normal *versus* 
#' $H_1:$ Os erros não seguem uma distribuição normal
#' 
shapiro.test(res_Stud)
#'
#' Como o valor-p é maior do que $\alpha=0,05$, não há evidências
#' para rejeitarmos $H_0$, considerando-se o nível de 5\% de 
#' significância. Portanto, a análise pode ser validada.  

#' ## Necessidade de transformação de dados
#' Como cultura geral, haja vista que as condições são atendidas

ggplot(dados, aes(x = fitted(modelo), y = res_Stud)) +
  geom_point() +
  geom_hline(yintercept=0, color = "red") +
  xlab("valores ajustados") +
  ylab("resíduos studentizados")

#' - Transformação Box-Cox
#' 
library(MASS)
with(dados,
     boxcox(y ~ trat,
            ylab="logaritmo da verossimilhança"))

#'## Obtendo o valor de lambda

box.tr<-boxcox(y+0.5 ~ trat, data=dados, plotit=T)
lambda <- box.tr$x[which(box.tr$y == max(box.tr$y))] 
lambda

#' Como $\lambda=1$ não há necessidade de transformação$
#'
#'

#' # Comparações múltiplas

#' ## Teste t com correção de Bonferroni
library(DescTools)
with(dados, pairwise.t.test(y, trat, p.adj = "bonferroni"));

#' ## Teste de Tukey
library(ExpDes.pt)
?dic
with(dados,
     dic(trat, y, quali = TRUE, mcomp = "tukey",
         sigF = 0.05, sigT = 0.05))

#' #Teste de Dunnett (sendo A o genótipo de referência)
library(multcomp)
summary(glht(modelo,
             linfct = mcp(trat =c("B-A == 0", 
                                  "C-A == 0",
                                  "D-A == 0"))))


#' #Outros testes aula dia 01/09 


#' ## Teste de Duncan
with(dados,
     dic(trat, y, quali = TRUE, mcomp = "duncan",
         sigF = 0.05, sigT = 0.05))






#' #Contrastes ortogonais
summary(glht(modelo,
             linfct = mcp(trat = c("A+B-C-D == 0"))))
summary(glht(modelo,
             linfct = mcp(trat = c("A-B == 0"))))
summary(glht(modelo,
             linfct = mcp(trat = c("C-D == 0"))))