#===============================================================================
# Título: Análise didática de experimento em parcela subdividida
# Curso : CET076 - Metodologia e Estatística Experimental
# Autor : José Cláudio Faria/UESC/DCET
# Data  : 06/04/2022 17:43:03
# Objetivos:
#===============================================================================
# O script foi adaptado e melhorado e partir do script feito pelo aluno
# de graduação ROGÉRIO SERWATKA ALONSO
#===============================================================================

#. Importante
# OBS: como o fator alocado na subparcela é quantitativo, o correto é analisar
# este experimento via análise quantitativa de experimentos, ou seja, regressão!

# Será mantido no material didático da disciplina apenas com finalidades didáticas!

#. Leitura dos dados
dados<-read.table('https://lec.pro.br/download/R/dados/eps.txt',
                  stringsAsFactors=TRUE,
                  header=T)

head(dados,
     13)
str(dados)

is.factor(dados$bl) 
is.factor(dados$adubo)
is.factor(dados$dose)
is.numeric(dados$y)

summary(dados)

#. ANOVA
av <-aov(y ~ bl + adubo * dose + Error(bl/adubo),
         data=dados)

summary(av) # Notar que a interação entre adubo e dose
            # é significativa

library(TukeyC)
cv(av)

#... Ftab. GL adubo 3; GL resíduo 9
(Ftab1 <- qf(5 / 100,
             3,
             9,
             lower.tail=FALSE))

#... Ftab. GL adubo:dose 6 ; GL resíduo 24
(Ftab2 <- qf(5 / 100,
             6,
             24,
             lower.tail=FALSE))

#... Ftab.GL dose 2; GL resíduo 24
(Ftab3 <- qf(5 / 100,
             2,
             24,
             lower.tail=FALSE))

#. GRÁFICO ADUBO
curve(expr=df(x, 3, 9),
      main="ADUBO - F calculado F(3, 9) gl",
      xlab="F",
      ylab="Densidade de probabilidades, f(F)",
      xlim=c(0, 10),
      n=1e3)

abline(v=Ftab1,
       col="red")

abline(h=0,
       lty=2)

xf <- seq(Ftab1,
          10,
          0.01)

ydf <- df(xf,
          3,
          9)

polygon(c(Ftab1, xf),
        c(0, ydf),
        col="red")

text(x=6,
     y=.1,
     paste("pf(>=",
           round(Ftab1, 2),
           ") =",
           round(pf(Ftab1,
                    3,
                    9,
                    lower.tail=F),
                 4)),
     cex=0.8,
     col="black")

#. GRÁFICO DOSE
curve(expr=df(x, 2, 24),
      main="DOSE F calculado F(2, 24) gl",
      xlab="Valor F",
      ylab="Densidade Probabilística F(f)",
      xlim=c(0, 10),
      n=1e3)

abline(v=Ftab3,
       col="red")

abline(h=0,
       lty=2)

xf <- seq(Ftab3,
          10,
          0.001)

ydf <- df(xf,
          2,
          24)

polygon(c(Ftab3, xf),
        c(0, ydf),
        col="red")

text(x=5,
     y=.1,
     paste("pf(>=",
           round(Ftab3, 2),
           ") =",
           round(pf(Ftab3,
                    2,
                    24,
                    lower.tail=F),
                 4)),
     cex=0.8,
     col="black")

#. GRÁFICO ADUBO:DOSE
curve(expr=df(x, 6, 24),
      main="ADUBO:DOSE - F calculado F(6, 24) gl",
      xlab="F",
      ylab="Densidade Probabilística F(f)",
      xlim=c(0, 10),
      n=1e3)

abline(v=Ftab2,
       col="red")

abline(h=0,
       lty=2)

xf <- seq(Ftab2,
          10,
          0.001)

ydf <- df(xf,
          6,
          24)

polygon(c(Ftab2, xf),
        c(0, ydf),
        col="red")

text(x=4,
     y=.1,
     paste("pf(>=",
           round(Ftab2, 2),
           ") =",
           round(pf(Ftab2,
                    6,
                    24,
                    lower.tail=F),
                 4)),
     cex=0.8,
     col="black")

#. INTERAÇÃO
par(mfrow = c(2, 1))

with(dados, 
     interaction.plot(adubo, 
                      dose, 
                      y, 
                      col = 1:4, 
                      lwd = 2, 
                      ylab = "Kg/ha"))


with(dados, 
     interaction.plot(dose, 
                      adubo, 
                      y, 
                      col = 1:4, 
                      lwd = 2, 
                      ylab = "Kg/ha"))

#. TUKEY
#.. Testando efeito principal das parcelas (adubo)
# (finalidade didática pois a interação é significativa: necessário desdobrar!)
adu <- TukeyC(av,
              which='adubo',
              error='bl:adubo',
              sig.level=.05)

summary(adu)

#.. Testando efeito principal da subparcela (dose)
# (finalidade didática pois a interação é significativa: necessário desdobrar!)
dos <- TukeyC(av,
              which='dose',
              error='Within',
              sig.level=0.05)

summary(dos)

#.. adubo/dose (dado adubo x, estudar dose)
#... adubo: nível_1
sali <- TukeyC(av,
               which='adubo:dose',
               error='Within',
               sig.level=0.05,
               fl1=1)  # 1: SALITRE C

summary(sali)

#... adubo: nível_2
sulf <- TukeyC(av,
               which='adubo:dose',
               error='Within',
               sig.level=0.05,
               fl1=2)  # 2: SULFATO

summary(sulf)

#... adubo: nível 3
urei <- TukeyC(av,
               which='adubo:dose',
               error='Within',
               sig.level=0.05,
               fl1=3)  # 3: URÉIA

summary(urei)

#... nível 4
caln <- TukeyC(av,
               which='adubo:dose',
               error='Within',
               sig.level=0.05,
               fl1=4)  # 4: CALNITRO

summary(caln)

#. CONTRASTES
# É bem trabalhoso no caso dos experimentos em parcelas subdividadas
# devido a mais de uma estimativa do resíduo!

# Matriz de contrastes dos adubos
                                      #A1  A2   A3   A4
mca <- rbind('A2 vs. (A1, A3, A4)' = c(-1,  3,  -1,  -1),
             'A3 vs. (A1, A4)'     = c(-1,  0,   2,  -1),
             'A1 vs. A4'           = c( 1,  0,   0,  -1))

# Matriz de contrastes das doses
                                  #D1  D2  D3
mcd <- rbind('D3 vs. (D1, D2)' = c(-1, -1,  2),
             'D1 vs. D2'       = c( 1, -1,  0))

library(gmodels) # Necessário para a função: make.contrasts
# Com a participação do prof. Ivan Bezerra Allaman
# foi possível fazer a decomposição da SQD desse experimento corretamente!

av3 <-aov(y ~ bl + adubo + adubo/dose + Error(bl/adubo),
          data=dados,
          contrasts=list('adubo'=make.contrasts(mca),
                         'dose'=make.contrasts(mcd)))

effects(av3) # Mensagem de erro -> método não aplicado a classe!
coef(av3)    # Dá para pegar a mesma informação com a função coef!

summary(av3)

summary(av3,
        split=list('adubo'=list('A2 vs. (A1, A3, A4)' = 1,
                                'A3 vs. (A1, A4)'     = 2,
                                'A1 vs. A4'           = 3),
                   'adubo:dose'=list('A1/dose*'          = c(1, 5),
                                     'A1/d3 vs. (d1, d2)'= 1,
                                     'A1/d1 vs. d2'      = 5,
                                     'A2/dose*'          = c(2, 6),
                                     'A2/d3 vs. (d1, d2)'= 2,
                                     'A2/d1 vs. d2'      = 6,
                                     'A3/dose*'          = c(3, 7),
                                     'A3/d3 vs. (d1, d2)'= 3,
                                     'A3/d1 vs. d2'      = 7,
                                     'A4/dose*'          = c(4, 8),
                                     'A4/d3 vs. (d1, d2)'= 4,
                                     'A4/d1 vs. d2'      = 8
                                     )))