#' ---
#' title: "Exemplo: Esquema fatorial de tratamentos"
#' author: "Renata Alcarde Sermarini"
#' date: ""
#' header-includes: \usepackage{here}
#' ---


#' # Croqui
#' 
#' Um possível croqui para o experimento fatorial $2\times 2$:
#' 
#' - Criando o plano de casualizaçãoseguindo o delineamento inteiramente casualizado 
#' com três repetições
#' 

library(agricolae)

Tratamentos <- c("I0:C0", "I0:C1",
                 "I1:C0", "I1:C1")

Plano.DIC <- design.crd(Tratamentos,
                        r = 3)

#'
#' -- Croqui (DIC)
#'

library(agricolaeplotr)
plot_design_crd(Plano.DIC,
                factor_name = "Tratamentos",
                labels = "Tratamentos",
                ncols = 3,
                nrows = 4)

#' - Criando o plano de casualizaçãoseguindo o delineamento casualizado em blocos com
#' três repetições
#' 

Plano.DBC <- design.rcbd(Tratamentos,
                         r = 3)

#'
#' -- Croqui (DBC)
#'

plot_rcdb(Plano.DBC,
          factor_name = "Tratamentos",
          labels = "Tratamentos")

#' - Criando o plano de casualizaçãoseguindo o delineamento quadrado latino (necessariamente com quatro repetições)
#' 


Plano.DQL <- design.lsd(Tratamentos)

#'
#' -- Croqui (DQL)
#'

plot_latin_square(Plano.DQL,
                  factor_name = "Tratamentos",
                  labels = "Tratamentos")


#' # Análise dos dados (exemplo)

#' Experimento fatorial 2x2 segundo o delineamento 
#' inteiramente casualizado com três repetições, 
#' para analisar o efeito da calagem e da irrigação 
#' sobre o peso de plantas.

#'
#' ## Entrada dos dados
#' 
rm(list = ls())
dados <- read.csv2("fat1.csv")
str(dados)

dados <- transform(dados,
                   Trat = as.factor(Trat),
                   I = as.factor(I),
                   C = as.factor(C))
str(dados)

#'
#' ## Gráficos
#' 
library(ggplot2)
ggplot(dados, 
       aes(x = Trat, 
           y = MS)) +
  geom_point()

with(dados,
     tapply(MS, I, mean))

plt.C <- ggplot(dados, 
                aes(x = I, 
                    y = MS, 
                    group = 1)) +
  geom_point(stat = 'summary', 
             fun = mean) +
  geom_line(stat = 'summary', 
            fun = mean) +
  labs(x = "Irrigação",
       y = "Massa seca")
print(plt.C)

with(dados,
     tapply(MS, C, mean))

plt.E <- ggplot(dados, 
                aes(x = C, 
                    y = MS, 
                    group = 1)) +
  geom_point(stat = 'summary', 
             fun = mean)+
  geom_line(stat = 'summary', 
            fun = mean)+
  labs(x = "Calagem",
       y = "Massa seca")

print(plt.E)

with(dados,
     tapply(MS, list(I, C), mean))

plt.int1 <- ggplot(dados, 
                   aes(x = I, 
                       y = MS, 
                       group = C, 
                       color = C)) +
  geom_point(stat = 'summary', 
             fun = mean)+
  geom_line(stat = 'summary', 
            fun = mean) +
  labs(x = "Irrigação",
       y = "Massa seca")
print(plt.int1)

library(gridExtra)
grid.arrange(plt.int1, plt.C, ncol=2)


plt.int2 <- ggplot(dados, 
                   aes(x = C, 
                       y = MS, 
                       group = I, 
                       color = I)) +
  geom_point(stat = 'summary', 
             fun = mean)+
  geom_line(stat = 'summary', 
            fun = mean)+
  labs(x = "Calagem",
       y = "Massa seca")
print(plt.int2)

grid.arrange(plt.int2, plt.E, ncol=2)

#'
#' ## Ajuste do modelo
#'  

modelo <- lm(MS ~ Trat, 
             data = dados)
#'
#'
#' ### Verificação das pressuposições da ANOVA
#'
res_Stud <- rstandard(modelo)
#'
#' - Normalidade dos erros
#' 
qqnorm(res_Stud)
qqline(res_Stud)
# library(hnp)
# hnp(modelo,
#     print.on = TRUE)
#'
#' Teste de normalidade:
#' 
#' $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)
#'
#' Considerando-se o nível de 5\% de significância não há evidências para rejeitarmos $H_0$. 
#' Logo, não se pode afirmar que os erros não seguem a distribuição normal.
#' 
#' - Homogeneidade de variâncias
#' 
ggplot(dados, 
       aes(x = Trat, 
           y = res_Stud)) +
  geom_point()

#'
#' Teste de homogeneidade de variâncias
#' 
#' $H_0$: há homogeneidade de variâncias *versus* $H_1$: não há homogeneidade de variâncias.
#'
library(lmtest)
bptest(modelo)
#'
#' Considerando-se o nível de 5\% de significância não há evidências para rejeitarmos $H_0$. 
#' Logo, não se pode afirmar que não há homogeneidade de variâncias.
#' - Relação entre média e variância
#' 
ggplot(dados, 
       aes(x = fitted(modelo), 
           y = res_Stud)) +
  geom_point()
#'
#' Box-Cox
#'
library(MASS)
boxcox(modelo)
#'
#' Não há necessidade de transformação dos dados, pois $\lambda = 1$ pertence ao intervalo de confiança.
#'
#' Como todas as pressuposições foram atendidas, temos o seguinte quadro da ANOVA.
#'
#'
anova(modelo)

modelo.fat <- lm(MS ~ I*C, 
                 data = dados)
anova(modelo.fat)

library(ExpDes.pt)
with(dados,
     fat2.dic(I,
              C,
              MS,
              quali = c(TRUE, TRUE),
              fac.names = c("Irrigacao", "Calagem")))